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CHAPTER I: INTRODUCTION 


The main objective of image data compression is to reduce the volume of data needed in 
order to represent an image. The major motivation for the study of compression, arises 
due to the restricted memory resources for data storage. The amount of data involved in 
Picture Archiving and Communication Systems (PACS) [1], is a typical example where 
data compression schemes can have a great economic impact. The advent of Internet and 
networking in general, has vastly increased the information transfer. 'The amount of data 
transmitted is limited by the available channel bandwidth. Data compression plays a 
major role not only in storage but also in transmission. Another useful application of data 
compression is the development of fast algorithms where the operations required to. 
implement an algorithm is reduced by working with the compressed data. 

Data compression can be achieved by exploiting the redundancy, which is a 
characteristic, related to predictability, correlation etc., in the data. The two main 
components of data compression are redundancy modeling and coding. In the first step 
(redundancy modeling), the data is transformed into a suitable domain where maximum 
information is packed in minimum number of samples. The objective of that 
transformation is to extract the redundancy present in the original data. Such a 
transformation will be data dependent, and will be based on the redundancy model. The 
second step (coding) actually stores the new data form with some predefined rules and 
conventions so as to obtain more compact data storage. 

In digital computers, data is represented as a string of binary digits, popularly known as 
bits. A bit can represent two distinct states 0 and 1. A string of k bits can be used to 
represent 2 k distinct states. Code is a string of bits assigned to each symbol of the data. 
Symbols in the data, which could be letters, signs, numbers and marks etc., are mapped to 
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codes and these codes are stored in digital computers. The original data can be recovered 
by using the relation between the symbols and codes. 

Let S = { Si, ,s m } be a set of m distinct symbols or states. A symbol is an atomic 

unit of information. The minimum number of bits used to code each symbol in S is the 

least positive integer k such that 2 k > m. If the data set D = {d j , ,dn} of length n, then it 

can be coded using nk bits. This kind of' coding is referred to as Fixed Length Coding 
(FLC), where the code length of each symbol in the data sequence is constant. 

Accurate model of the data redundancy, helps to characterize data sets based on the 
considered redundancy. Compression algorithms based on the redundancy models are 
very effective for the class of data characterized by certain measure of redundancy. The 
objective of this work is centered on the compression of image data, particularly for the 
images arising from Magnetic Resonance Imaging (MRI). Medical imaging in general 
aims at producing images of sections of a physical object for diagnostic purposes. MRI is 
one of the most prominent medical imaging modality. MRI images have better soft tissue 
contrast than the images produced using Computerized Tomography (Hermen, [27]). 
Another important area where data compression schemes are highly applicable is the 
Teleradiology (Gitlin, [22]), or the practice of radiology at a distance. It involves 
transmission of medical images over a distance. Not much loss of information due to 
compression and decompression can be tolerated in most medical applications. These 
stringent requirements are the motivation for designing sophisticated compression 
algorithms for medical images. 

1.1 Digital Image Model: 

An image in general is an analog intensity distribution function which can be viewed as a 
mapping from [0,2 tc] x [0,27t] to [0,1]. This continuous intensity distribution function is 
digitized for storing the image in a digital computer. Let f be the analog intensity 
distribution function given by. 
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Consider for 0 < i<N, 0< j<N, 


Xij(x,y)=- 


N 2 

0 


f : [0,2 tc] 2 -> [0,1] 


■r i . i+1 

if — < x < 

N N 

else 


— < y< 
N 


1+1 

N 


Let, 

P ij (0= JJf (x» y)X fl (x, y)dxdy , where I = [0,2rc] 2 

I 

and f « F(x, y) = -h 2 P y (f )X.. (x , y) . (1.1.1) 

i,j 

Let F be the digital image. Py ‘s are the pixel values to be stored for representing the 
digital image. 


Depth and Resolution of a digital image: 

Resolution of a digital image is the discretization of the domain [0,2tc] x [0,2tc]. If N x N 
is the discretization of the domain where N = 2 m for some m e Z+, then the resolution of 
the digital image is m. Cases of the most practical interest are m = 7,8,9,10. 

Depth is the discretization of the range [0,1] and it is also the number of distinct gray 
levels the image can take. It is expressed in bits per pixel (bpp). An image is of b-depth 
if the range [0,1] is discretized into 2 b parts. 

Resolution, and depth are the two quantities that determine the storage requirements. 
Throughout this thesis, by data we mean image data. The depth of the image would 
determine the number of states or symbols an image can take. Some useful relationship 
between the depth and resolution of a digital image is obtained in the fourth chapter of 
the present thesis. 
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1.2 Variable Length Coding (VLC): 

In 1952, Huffman [30] gave a coding scheme where the code length is not constant for 
each symbol. It is based on the idea that the frequently occurring symbols can be coded 
with codes of smaller length as compared with the less frequently occurring symbols. 
The success of this coding in terms of compression is based on the hypothesis that all 
symbols in a data sequence are not equally likely. It is necessary for such codes to be 
uniquely decodab le. Before going into the details of VLC, let us introduce some standard 
definitions, which will be used frequently. 

Frequency distribution of the symbol set S consisting of m symbols, in the data D of 

length n, is the collection of positive numbers, Fr = {fi, ,f m } one for each state, defined 

by fj = # { dk € D : dk = Sj } , where # denotes the cardinality of the set. It follows that 

m 

fi = n . In 1940 c s, Shanon [68] defined the zero order entropy of a data sequence D to 

i=l 

be 

m 

e = e(D) = -^p i log 2 p i , (1.2.1) 

i=l 

where pj is the probability of occurrence of i th symbol in the data sequence D. Let 

C = {ci, ,c m } be a code book. It follows that the average code length C, to represent 

the data sequence D is e c (D) given by, 

m . . 

e c (D) = ypj|Cj|, where \c i | stands for the length of the code Cj. 

i=l 

Theorem 1.2.1 (Shanon Entropy Theorem): 

Given a data sequence of symbols, for any code book given to that symbol set, the 
average code length using the code book to represent the data sequence is always greater 
than or equal to the zero order entropy of the data sequence. 

m m 

i.e. e c (D) = £Pi| c i| > e(D) = -£p ; log 2Pi . 

i=l i=l 
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We state the definition of Markov process, before describing certain algorithms. These 
algorithms try to achieve minimum code length prescribed by Shanon [67]. 

Let {d n } be a data sequence generated in some process. It follows a k-th order Markov 
model (see, Sayood, [65]) if the knowledge of past k-symbols is equivalent to the 
knowledge of the past history of the process. 

= f(dj 

The values taken by (d n _, , ••• , d n _ k ) are referred to as the states of the process. For an 

image of depth b, the number of states the pixel values can take is 2 b . In the sequel we 
will use 2 b as the number of elements in the symbol set S. 

1.2.1 Huffman Coding: 

Huffman [30] proposed a method for reducing the code length for data sets in which each 
symbol occurs with a different probability and has given an algorithm for constructing 
variable length codes. We present the algorithm here. 

Consider, a zero order Markov process, given by a symbol set, S = js, ,• , s 2 „ j . In the 

case of images, elements of the symbol set S will contain gray level values. Here S 
denotes set of states an image of b bpp depth can take. And Fr is the frequency set given 

by Fr = jf,, ,f 2 „ j , where f, denotes the number of occurrences of the symbol s*. 

Given the data set this quantity is calculable. This is a redundant information present 
along with the data. „ ^ 

Coding: 

Step 1 : We pick two symbols from the set S which have frequencies lower than the other 
symbols in the given data set D. Let the chosen symbols be denoted by c and d. 

Let cr = a 2 " : {l,2,- • ■ ,2 b } -» {l,2,- • • ,2 b } denote a permutation such that 
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c = s a(2 b -i) ’d = s 0 ( 2 ‘> ) - Then we define a new Markov source consisting of the set 
of symbols, S 2 = js c(1) ,---,s c(2b _ 2) ,t 2 ~‘j and the set of frequency distribution, 

Fr 2 1 = |f<,(i ) >‘"»f 0(2 k_ 2) >f c(2 »_, ) +f j(2 b ) |-The 2 b state problem has reduced to a 

2 b ’’ states problem of finding a code book C 2 "' 1 = jc 2 '’~ 1 ,-‘ 4 ,c 2 b _' ] j for S 2 ”- 1 . In 
summary, the first step can be described as 
{S,Fr,c} => cr b => {s^Fr 2 *- 1 ,^- 1 } . 

Step 2: Step 1 is applied to the 2 b - 1 symbol Markov source, |s 2b " 1 ,Fr 2b_I ,C 2k_I j . So 

that the 2 b -l symbol problem is reduced to 2 b -2 symbol problem. 

Step 3: Step 2 is repeated and after 2 b -l steps, we obtain {S 2 ,Fr 2 ,C 2 } a data source of a 
symbol set S 2 , that has only two symbols in it. 

Step 4: To the symbol set S 2 , which has two elements, we assign a code of 0 to one 

symbol and 1 to the other symbol, i.e. C 2 = {c 2 = 0,c 2 = l] . 

Step 5: We build the codes C 3 = |c 3 ,c 2 , causing, c 3 3(i) =c 2 , c 3 3(2) = c 2 0, c 3 3(3) = c 2 l. 
Here, c 2 0 and c 2 l denote the string c 2 concatenated with 0 and 1 respectively. 
Finally, codes C = |c ] ,c 2 ,---,c 2b _ i | are obtained using the relations c o „ (j) =c 2 _1 
for i = l,2,....,2 b -2, C CT „ (2b _ l} =c 2 b "|0 and c a „ (2b) =c 2 b “|l. 

Either the code book or the frequency set Fr is stored along with the coded data. In case 
the code book is not stored, it can be constructed by using Fr. 

Decoding: 

Step 1: We read one bit of the coded data and store it as a string bi. 

Step 2: String bi is matched with the code book, if a matching occurs bi is decoded, else 
step 3 is executed. 

Step 3: Another bit is read and concatenated with the string bi. Step 2 is repeated until 
all the bits in the coded data is exhausted. 
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Length of Huffman codes: 

Let 1 h be the average code length in a Huffman coding for a source S. Let, E(S) be the 
entropy of the source S, given by equation (1.2.1). Gallager [20] has shown that the 
average length of Huffman code satisfies the following inequality, 

2(log 2 e) ~ 

e j 

where pi is the probability of most likely symbol in the data sequence. There are several 
modifications to the Huffman coding scheme available in the literature (Kunth [43], 
Sayood [65]). A combination of pixels considered to constitute one symbol is one such 
modification to the original Huffman scheme. The average bit rate using modified 
Huffman coding, approach the true entropy of the image, if combinations of pixels are 
taken to constitute single symbol. In Huffman coding, we require integral number of bits 
to represent each symbol. Even if the histogram of an image is highly peaked, this 
coding would require at least one bit per pixel. Apart from lossless image compression, 
Huffman coding has applications in text and audio compression. 

1.2.2 Arithmetic Coding: 

Arithmetic coding is another popular method for data compression, and it does not suffer 
from the limitation that each code word has to be one bit long. There have been many 
developments in the practical arithmetic coding strategies by Langdon and Rissanen [44], 
Here we present the basic algorithm for arithmetic coding. 

Let S = { 1,2, ,2 b } be the source set of 2 b symbols, typically representing an image of 

depth b. Let 

f c :{l,2,--,2\2 b+I } [0,1] 

be the cumulative distribution function (cdf) of the probabilities P = jpj,---,p 2b j of 
occurrences of the 2 b symbols in a digital image. 

f c( i ) = EPi> for i = l,2,---,2. b . 
j<i 

Let us assume, f c (0) = Oand f c (2 b +1) = 0. 


E(S)<l H <E(S) + p I+ log 2 
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Let D = {d, , • , d m | be a data sequence of length m. We need to find the code which 

is a real number here. The real number obtained in this process is stored in sufficient 
precision so that the data sequence can be uniquely decoded. The following steps are 
used, 


Coding: 

Step 1: Set x 0 = 0, y 0 =1. Fork =1,2, ,m+l, the values of Xk and yk are 

calculated using the relation, 

X k = X k-1 + (Yk-! - X k-1 ) f c ( d k " 1) 

y k = x k-. + (y k -i - x k-i) f c ( d k) 


( 1 . 2 . 2 ) 


Step 2: C = 


X m + y n 


, where C is the code for the input data sequence. 


Decoding: 

Q — x ■ 

Step 1: We set x 0 = 0 and y 0 = 1. For each i, we find C’ = Lli — 

y i - 1 ~ X i-1 

Step 2: We find the value of dk for which C* e [f r (dk-l), f(dk)). The values of Xk and yk 
are updated using the relation (1.2.2). 

Step 3: Step 2 is repeated until the entire data is decoded. 

Length of Arithmetic codes: 

The average length of an arithmetic code U of a data sequence of length m is given by, 

E(SJ < 1 A < E(S m ) + 

where E(S m ) is the entropy of sequence of m symbols. By increasing the length of the 
sequence, a rate as close to the entropy as we want can be guaranteed. Arithmetic coding 
method is used in Joint Bi-level Experts Group (JBIG) standard (Langdon [37], Sayood 
[65], Hunter and Robinson [31]) for encoding binary (two color) images. 
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1.2.3 Dictionary Based Compression: 

It is one of the approaches to encode data sources in which certain patterns are frequently 
occurring. The encoding method is to keep a dictionary of frequently occurring patterns. 
When these patterns appear in the source data they are encoded with a reference to the 
dictionary. The principle of this kind of compression is that it encodes variable length 
strings of symbols as a single token, which are some reference to the dictionary. The 
objective is to replace certain phrases with tokens. If the number of bits in the token is 
less than the number of bits in the phrase, compression will occur. There are two 
approaches for building the dictionary. Static dictionary, where the dictionary is 
completely defined when encoding begins. This approach would require sufficient prior 
knowledge about the data. On the other hand, Adaptive dictionary start out with a default 
baseline dictionary. As the coding proceeds, the algorithm adds new phrases to be used 
later as encoded tokens. 

Certain aspects of statistical modeling, related to entropy, character and word 
frequencies, were the main areas of research in data compression up to 1977. In the years 
1977 and 1978, J.Ziv and A.Lempel [77]-[78], initiated works in dictionary based 
compression. They developed two algorithms LZ77 and LZ78, and these algorithms 
triggered research in dictionary based compression. In principle, these are two different 
techniques. 

LZ77 approach: 

The input data sequence is scanned by a sliding window. The window consists of two 
parts, a search buffer that contains the portion of the recently encoded sequence, and a 
look-ahead buffer that contains next portion of the sequence to be encoded. The sizes of 
buffers are implementation dependent and they are generally large. The algorithm moves 
a search pointer back through the search buffer until it encounters a match. The number 
of symbols passed by the pointer to reach the match is called offset. The algorithm then 
examines the symbols following the symbol at the pointer location to see if they match 
consecutive symbols in the look ahead buffer. Match length 1, is the number of 
consecutive symbols in the search buffer that matches consecutive symbols in the look 
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ahead buffer. After finding the longest match for the given window size, the algorithm 
codes it with a triplet { o , 1, c} where o is the offset, 1 is the match length and c the code 
word corresponding to the symbol in the look ahead buffer that follows the match. 

LZ78 approach: 

In LZ78, instead of using fixed length phrases from a window into the data sequence, it 
builds phrases up one symbol at a time, adding a new symbol to an existing phrase when 
a match occurs. It constructs a dictionary while both coding and decoding. The 
construction at both ends is done in an identical manner. The inputs are coded as a 
double {i, c} where i corresponds to the index to the dictionary entry that was the longest 
match to the input, and c being the code for the character in the input following the 
matched portion of the input. In case of no match the double becomes a new entry in the 
dictionary. Each new entry into the dictionary is one symbol concatenated with an 
existing dictionary. 

The dictionary building is more implementation based and several variations are possible. 
Notable among the modifications is the. one given by Welch [74], and LZ is modified to 
LZW algorithm. 

Run Length Coding: 

It is a simple technique to compress large runs of identical symbols in a data sequence, 
(Capton [11]). RLE encodes a run of symbols as a symbol and a count. It is often used 
for its easy implementation. 

The pixel values in most practical images tend to be well spread out over their entire 
range. The peaks in histogram may not be found in images for statistical compression to 
succeed. 
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1.3 Lossy Compression: 

The major advantage of the above methods is that, after decompression we get back the 
original image without any distortion. These types of algorithms are in general referred 
to lossless compression. These methods do not give rise to high compression when 
applied to image data. Images based on scanned photographs dp not posses the 
characteristics to create multiple occurrences of the same phrase. Out of a sequence of 
twenty-five pixels say, one or two may vary by a single step from the scans before and 
after and they often reduce the performance of dictionary or statistically based 
compression. However, these differences are small enough, to the extent that they are 
either undetectable or meaningless to human eye. We see that the pixel values 
necessarily have to follow similar pattern for the dictionary methods to work better. The 
minute variations which are undetectable to the human eye limits the effectiveness of 
these compression strategies. Certain mechanism which will enhance the performance of 
the above said methods are needed. 

In order to enhance the performance of the statistical or dictionary methods, models 
based on human visual systems (HVS), (Jayant, etal. [36]) were brought into the scene. 
In these models, details which are not detectable by human eye are considered as 
redundant. Consequently methods which introduce certain distortion between the 
original image and the decompressed image are developed. The algorithms are modeled 
in such a way that the distortion cannot be perceived by human eye. These algorithms 
are based on the philosophy that images are meant for visual inspection by humans, and 
they cannot perceive certain very fine variations in the pixel values, so it is sufficient to 
store images up to details which are perceivable by human visual system. These kinds of 
compressions are referred to as lossy compressions. However, the term ‘loss’ in the 
lossy compression depends on the application. Medical images helps radiologist, in 
diagnosis, it often imposes stringent conditions on the amount of permitted distortion for 
the sake of compression. Consequently lossy compression schemes for medical images 
would be required to perform in a ‘nearly lossless’ maimer. The general principle of this 
type of compression is predictive or transform coding followed by quantization. The 
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quantized transform coefficients are retained, from which the decompressed image could 
be reconstructed at the time of need. 

1.3.1 Predictive Coding: 

Methods which predict the current pixel (the one to be encoded) reasonably well by using 
the information about the previous pixel values fall in this category. It is based on the 
statistical dependency of one sample element and the next. Let {x n } represents a data 
sequence in which the first k-1 elements are encoded and we are encoding the k th 
element. Let x* be a quantity predicted for the data point x k based on the previous k-1 
elements, considering x k as an estimate of x kj the prediction error ek is defined as, 

e k = x k - -x*. 

It is sufficient to quantize e k instead of x k the actual data point. If e* is the quantized 
value of e k , then an approximation the original k-th data point x k could be obtained 
using, 

x k = x* + e* . 

Common data compression that uses predictive quantization is called Diffemtial Pulse 
Code Modulation (DPCM) (see, Jain, [35]). The error or loss in the process is given by, 

x k - x k = e k - e k . 

1.3.2 Transform Coding: 

Image data is transformed into a suitable transform domain by using an invertible 
transformation. Let the digital image be, (cf: 1.1.1) 

F = Z P ij( f ) X ij( x ’y)- 

U 

Pij(f) are stored as pixel values of the image f. A representation of f in the transform 
domain is obtained as f , given by 

f = 0-3.1) 

ij 

where C ;j (f) are the transform coefficients of f and {<j)jj }is an orthonormal basis of the 
transformation. In transform coding, a discrete transform based on is applied over 
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the pixel data to get the transform coefficients. Instead of storing the transform 
coefficients Cjj(f), quantized coefficients C*j(f)are stored. C ^(f) ’s are obtained by 

quantizing C s j (f) ’s. Some of the useful candidates for the transform are listed here. 
Karhunen-Loeve Transform (KLT): 

It is the optimal transform coder (see Donoho, et al. [19]), in the sense that it minimizes 
the mean square error of the decompressed data for a given number of total bits. 
Consider a random vector, 

X = [X,, ,X n f with zero mean. The vector X is linearly transformed by a matrix A 

to produce a vector Y = [Yi, ,Y n ] T . The auto covariance matrices R x and R y are 

defined as, 

R x = E (XX T ) and R y = E(YY T ). 

The transformation matrix A is chosen such that the components of the resulting vector Y 
are mutually uncorrelated. After obtaining Y, it is subjected to component wise 
quantization to obtain Y an approximation to Y. The vector Y , is linearly transformed 
by a matrix B to get X . Matrices A and B are obtained so that the quantity, 

The matrix R x is symmetric and positive semi definite, there exists a full set of eigen 
vectors with non negative eigen values. The KL transform matrix A for the data 
sequence X is constructed by keeping the orthonormalised eigen vectors of the auto 
covariance matrix R x . The rows of R x are ordered in terms of decreasing eigen values, 
that is R X A=AA where A = diag(A.i,....,X n ). Transforming X with A will diagnolise R y . 

R y = E[ A* XX t A] = A* E(XX t ) A = A* R x A = A. 

KLT satisfies the best linear approximation property in the 1 2 sense. When the number of 
transform coefficients to be retained are fixed, then KLT is the optimal transform. 
However the algorithm complexity involved in KLT, restricts its practical usage. It 
requires N 2 operations per vector of length N. If we view the image as a N X N matrix, 
comprising of N row vectors or N column vectors, then KLT would require N 3 operations 
on these N vectors. Also KLT is dependent on the data sequence to be compressed. Due 


E (X(i) - X(i)) 2 


,i = l 


is minimized. 



Chapter 1 


14 


to these reasons, fast fixed transforms, which work at N 2 logN complexity, are practically 
used for compression. Discrete Fourier transforms, Discrete Cosine transform, Discrete 
Sine transform, Hadamard transform etc., are being substituted for KLT, because of their 
implementation efficiency. The following transforms are commonly used (Jain, [34]). 


Discrete Fourier Transform: 

The DFT forward and inverse for an image is given by, 

1 N-l N-l i(mx+ny)27t 

DFT(m,n)=^££Pix(x,y)e" N 

-IN X _Q y_0 

i N-l N-l i(mx+ny)27C 

Pix(m,n) = — J]J]DFT(x,y)e N . (1.3.2) 

-IN x =0 y=0 

Although the Fourier transform produces complex transform coefficients in general, for 
real input sequences, it is sufficient to store half the data and the rest can be generated by 
using the symmetry properties. 


Discrete Cosine Transform: 


N-l N-l 


DCT(i, j) = -C(i)C(j)XZPix(x,y)cos 


x=0 y=0 


pix(i, j) = ^XE G ( x ) c (y) DCT ( x ’y) cos | 


(2x + l)rri 
2N . 

(2x + l)7ii 
2N 


cos 


cos 


(2y + l)7ij 
2N 

(2y + l)rtj 
2N 


where C(x) = ( 



if x = 0 
if x > 0 


(1.3.3) 

(1.3.4) 


DCT performed on an N x N square matrix of pixel values yields an N x N square matrix 
of frequency coefficients. DCT can be actually computed like FFT and can be done in 
0(N 2 logN) complexity. 


Discrete Sine Transform: 

The sine transform and inverse sine transform pair for each row of the image is given 
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DST(j,k) = 


2 N “* N_1 

~2]^pix(m, n ) sin 

* m=0 n=0 


N + 


rt(j + l)(m + 1) 7r(k + l)(n + 1) 


N + l 


Pix(m, n) = — •£ £ DST(j.k) sin 

VN + lJJfcJ N + l 


sin 


N + l 

0<k<N-l,0<j<N-l 

7i(k + l)(n + 1) 

N + l 

0<m<N-l,0<n<N-l 


(1.3.5) 


Hadamard Transform: 

HT(j, k)= -L |]|]pix(m,n)(-l) b(k ’ m)+b(k ' n) 0 < j < N - 1,0 < k < N -1 

v N m=0n=0 ^ ^ ^ 

Pix(m, n) = 1 gg H T(j,k)(-l) b(lc ’ m)+b(k ’ n) 0<m<N-l,0<n<N-l 

VN j=o k =0 

where, 

n - 1 

b(k, r) = ^ kjTj with k i5 r, = 0,1 and n = log 2 N. {kj} and {r,} are the binary 

i = 0 ' 

representations of k and r, respectively, that is, 
k = ko+2ki+...+2 n ’k n -i and r = ro+2ri+...+2 n *r n .i. 

The above mentioned transforms falls in the category of fast transforms. These 
transforms can be calculated via FFT, and their computational complexity is 0(N 2 logN) 
for N x N images. These transforms posses good energy compaction property for highly 
correlated images. DCT is a most widely used transform for image compression. 

1.3.3 Quantization: 

After the application of transform, the coefficients are subjected to quantization which 
actually causes the compression. Quantization, (Gray, R.M., etal., [23]) or round off is a 
technique originated from analog to digital conversion of signals. It amounts to 
partitioning the range into small sub intervals (subsets) and choosing a representative 
value for each sub interval. Any data point that occurs from the subinterval is quantized 
by the representative value. More precisely, a quatizer is a set of intervals or cells, S = { 
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Sj : i e Z + }. Together with a set of representative values C = {y; : i e Z + } so that the 
overall quantizer q is defined by 

q(x) = Zy ; X Si (x), (1.3.7) 

i 

where X s . (x) is the indicator function for the cell s,. S is the partition of the original 

range. Uniform quantizer is one in which the levels are equi spaced and the representative 
value is the mid point. 

The goodness of reproduction measures the quality of a quantizer. If x is the input signal 
and x is the output signal after quantization, the error due to quantization in 1 2 sense is 
given as 

q-l 2 = E[(x- x) 2 ). 

The quantizer that minimizes the 1 2 error for a given number of quantization levels is the 
Lloyd-Max quantizer (Jain, [34]). The concept of scalar quantization, where the 
quantization is done component wise has been extended to vector quantization (Gersho 
and Gray [23]), it is just like scalar quantization except that all components of a vector, 
say k successive source samples are quantized simultaneously. They are characterized by 
a k-dimensional partition, a k-dimensional code book (containing k-dimensional points, 
reconstruction code words or code vectors) and an assignment of binary code words to 
the cells of partition. 

1.4 Standard Compression Methods: 

In the present section, we give a brief introduction to the standard compression methods. 
These methods are compared with our strategies in the thesis. 

1.4.1 JPEG Standard: 

Constructive Committee on International Telephone and Telegraph (CCITT) and 
International Standards Organization (ISO), in 1980’s constituted a committee - Joint 
Photographic Experts Group (JPEG), with an idea to develop a standard for image 
compression. The committee came up with JPEG standard, which has been recognized 
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and used as a major portable graphical file format (Wallace, [73]). The JPEG committee 
designed a standard for compressing either full color or gray scale still images of natural, 
real world scenes. JPEG is designed to exploit the known limitations of the human eye, it 
falls in the category of lossy image compression. However, the standard allows the user 
to choose the amount of permitted loss. The degree of loss could be adjusted using a 
compression parameter called quality factor. The JPEG specification is broadly divided 
into lossless and lossy encoding modes. 

"The lossless compression uses the predictive or adaptive model described earlier, with a 
Huffman code output. It does not allow any distortion in the decompressed image. The 
compression efficiency in terms of the size of the compressed image is not satisfactory 
for many applications. On the other hand, the JPEG lossy compression algorithm has 
good compression efficiency. JPEG has three major parts - transform coding using 
Discrete Cosine Transform (DCT) (1.3.3)-(1.3.4), quantization and loss less compression 
of quantized DCT coefficients. 

DCT applied to the pixel data prepares the data for quantization. The amount of 
quantization will depend on the required quality of the output. Coefficients which are 
sufficiently close to zero are truncated to zero. The quantization step makes the DCT 
coefficients matrix sparse. Run length coding (RLE) is used to absorb the zero runs of 
the coefficient matrix. In order to increase the length of runs, the coefficient matrix is 
first rearranged to form a zig-zag sequence. Instead, of compressing the coefficients in a 
row major fashion, the zig-zag scanning moves along diagonal paths which progressively 
moves from low frequency to high frequency terms. Since the probability of a high 
frequency term rounding off to zero in the qunatization step is more, the coefficients are 
arranged in a frequency order so that the indexing will be minimal and easier when 
storing the quantized DCT coefficients. 

Finally an entropy encoder (either Huffman coding or Arithmetic coding) depending on 
the choice of the implementation is used. To enhance the speed of operations block 
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coding is done. In block coding, the whole image is partitioned into 8x8 blocks and the 
DCT, quantization and coding are applied for each block independently. 

Although JPEG works well for images of natural and realistic scenes. It does not work 
well for line drawings, cartoons and sharp edges in images tend to come out blurred. 
Human visual system is more sensitive to brightness variations than to hue variations. 
JPEG can compress hue data more heavily than brightness (gray-scale) data. Further 
JPEG accumulates losses, in other words recompressing a decompressed image loses 
more information. 

Another generation of image coders mostly based on wavelet decompositions, other fast 
transform based decompositions, different quantization techniques and entropy coding 
are being considered for the next standard JPEG-2000 (see Donoho, et al. [19]). 

1.4.2 Graphics Interchange Format (GIF): 

It is another popular commercial compression standard (Sayood, [65]). It uses LZW 
compression strategy (Welch, [74]) and certain variations of that. It is suitable for images 
that contain large regions of constant value. GIF does significantly well on images, 
which have few distinct levels. GIF is lossless for images having 256 distinct states (gray 
levels), which is the case of most practical interest in medical imaging. 

1.4.3 Fractal Compression: 

Objects which have self similar structures are called fractals. The main idea behind the 
fractal compression is to use the similarities of certain regions of object with itself. A 
mapping f from a set to itself is said to be contractive if d(f(x),f(y)) < d(x,y) where d is 
distance measure. It is well known that a contractive mapping has a fixed point (Rudin, 
[64]). Suppose that for a given image I, there exists a function f such that f(I) = I. Then f 
can be used for representation of the image I. Fractal compression amounts to finding a 
function f for the given image I, so that I is a fixed point of f. Arnaud Jacquin [33] in 
1989 came up with an idea of partitioning the image I into sub images Ik, and for each I k a 
suitable function f k is obtained by using a method called Iterated Function System (IFS) 
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(Y. Fischer [76]). These types of compression strategies are of recent origin and a lot of 
research needs to be done based on these ideas. 

1.4.4 Wavelet based Compression: 

In 1992,. R.A. DeVore, B.Jawerth and B.J.Lucier [17] have proposed a compression 
algorithm based on wavelet decomposition. They associate to each set of pixel values a 
function f given by, 

f = Z , (1.4.1) 

ka0,j=(j,,j 2 ) 

where c Jjk ’s are the wavelet coefficients of f. (]>j, k ’s are generated from a single function (j) 
by their dyadic dilates and translates given by, 

Mx) = <K2 k x- j)- (1-4.2) 

After obtaining a wavelet-based decomposition of the form (1.4.1), their compression 
algorithm quantizes the transform coefficients to obtain quantized coefficients c j k and 

the compressed function takes the form, 

f= Ic*. (1-4.3) 

kSO, j=(j, , j 2 ) 

Their theory relates the smoothness of images in certain smoothness classes called Besov 
Spaces [16], which are defined as follows: 

Definition (1.4.4) (Besov Spaces): 

Fix a > 0 and an integer k > 0. Define the k-th forward difference with a prameter h e 
* + by 

A° h f(x):=f(x) 

and for j = 1, ,k A j h f(x):= A J k ' f(x + h) - A J k ' f(x). 

For 1 < p < oo, the k-th modulus of smoothness in Lp(I) is defined as 
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® k (f,t):==sup(|A k h f(x)|). 

Nst 

where I r h = {x el: x+rh el}. A k h f(x) is defined iff x e I rh . The Besov space, 
B“ ,k (L p (I)) norm is defined as 


'Bj-‘(L p (l)) 




q 


for 1 < p < oo, 0<q<oo and for q = oo, 

II f ilB- k (L p (I)) := ll f llp+SUp[r a CO k (f,t) p ] 

f e B“- k (L p (D) if ||f|| B .. k(Lp(I)) 


The error function in Lp(I) is given by, 


a N (f) = inf f-f 

^ f hasSN coefficients ”L p (I) 


(1.4.5) 


For certain choices of cp, DeVore, Jawerth and Popov [18] have shown that. 


a N (f)n = 0(N"“ ) 


<=> 


f e B“’ (L (I)) 


(1.4.6) 


all 
where — + — = — 
2 p q 


They have proposed a compression algorithm based on the wavelet decomposition of f. 
The quantization strategy to obtain c j>k ’s is that they satisfy. 


( c « - £ 


N 


-l/q 


(1.4.7) 


Here N, the quantization parameter, is a chosen positive integer, and it determines the 
amount of compression. The number N, of nonzero coefficients c Jjk ’s is shown to 

satisfy, 

NSC.NIC.,,,,,, . 


(1.4.8) 
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1.5 MRI image compression: 

Images arising from various sources form different classes. For example images of 
natural scenes, machine parts, buildings, person’s etc., fall in different categories of 
images and the corresponding redundancy models are different. In the present thesis, we 
deal with the images arising from MRI. Apart from comparing a number of image 
compression methods, we propose newer methods for compression, which are well suited 
for MRI images. We study the reliability of the proposed compression methods. The 
performances of these methods are evaluated theoretically and experimentally. The goal 
of compression is to represent MRI images with the smallest possible number of bits, 
thereby minimizing storage requirements for archiving and speeding transmission for the 
sake of teleradiology (Gitlin, [22]). A compression system typically consists of one or 
more of the following units, 

(i) Approximation of the underlying intensity distribution function, by a suitable 
approximation process. 

(ii) Quantization of the analog or high rate digital pixels, into a relatively small 
number of bits. This conversion is lossy. 

(iii) Lossless compression achieved by invertible methods such as Huffman coding, 
Arithmetic coding, LZ coding (cf: section 1.2). One of the lossless compression 
methods is applied over the quantized data for obtaining net compression. 

The above units are studied for MRI images, in which not much loss of information can 
be tolerated in compression/decompression for most medical applications (Cosman, etal., 
[13]). Most of the existing methods produce compressed images, which are to be used, 
for visual inspection, and they use the limitations of Human Visual System (HVS). 
However, for MRI images, there are certain other quantities that can be computed for the 
whole image or for some sub regions, depending on the need. These quantitative 
informations along with the images are required in diagnostic studies, especially in 
cancer tissue characterization and to study the effectiveness of the treatment. Images 
produced using our compression strategy can be used for one or more of the following 
applications: 
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(cl) Estimation of tissue parameters, p, Ti and T 2 (studied more precisely in the next 
section) which are used for certain diagnostic studies. 

(c2) Volume or area estimation of geometric structure typically cancerous portions seen 
in MRI images. 

(c3) Visual inspection by experts, for diagnostic purposes. 

We note that the computations involved in (i) are highly sensitive to the input data. They 
require large set of input data and hence compression will be useful. Certain MRI 
preliminaries are presented in the following sub section. 

1.5.1 Magnetic Resonance Imaging - Basics: 

Magnetic Resonance Imaging (MRI) is a technique for providing a picture of the interior 
of a physical object. Magnetic resonance (MR) is based on the interaction between a 
system of atomic nuclei and an external magnetic field. In MRI, the system is perturbed 
from the equilibrium and the system relaxation to the equilibrium is observed. An 
induced electric signal is measured during this relaxation process. This signal depends on 
four quantities p, Ti, T2 and P (a parameter set), which measure the nuclear spin density, 
the interactions of the nuclei with their surrounding molecular environment, those 
between close nuclei, and certain pulse sequence parameters respectively. Based on a 
suitable set of measured electric signals, a function of p, Ti, T 2 and P is computed, which 
is called MR image. Basic references describing MRI are available in (Healy, et ah, [26], 
Hinshaw, et al., [29], Kumar, A., et ah, [42], Mansfield, et ah, [51], [52], Sebastiani et 
ah, [66], Stark et ah, [69]. Here we briefly present certain MRI preliminaries, which are ■ 
used in the present thesis. 

The hydrogen nuclei, often referred to as spins, can be imagined as tiny bar magnets. In 
the absence of an external magnetic field, spins assume random orientations, and the net 
magnetization will be zero. When exposed to a static magnetic field Ho, the randomly 
oriented nuclei lineup with the magnetic field. The spin vectors experience a torque or 
couple when subjected to the magnetic field, analogous to a spinning top in the earth’s 
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gravitational field. As a result they precess around the axis of the magnetic field at a rate 
given by the Larmour relationship (see Stark et ah, [69]), 

co yH 0 


f = 


(1.5.1) 


2tc 2 % 

where f is the resonance frequency in Hertz, co is the angular frequency in radians per 
second, y is the gyromagnetic ratio, which is the characteristic of every isotope and Ho is 
the static magnetic field. Let the vector sum of the individual spins be approximated by 
the net magnetization vector of the sample, M. The net magnetization M, is the actual 
observable property of the sample. 


In a typical MR experiment, the direction in which, the net magnetization vector M is 
pointing is disturbed. Then it relaxes back to its equilibrium position. The magnetic field 
Ho = Hok is applied in the z-direction, k is the unit vector in the z direction. Before 
excitation, M is directed along the z-axis and therefore has no component in the xy plane. 
The excitation process is tipping M away from the equilibrium position along z. This is 
achieved by temporarily applying a second, much weaker magnetic field Hi, 
perpendicular to Ho. Also Hi rotates about Ho at the Larmor frequency. Ho is normally in 
the range 0.1 - 5.0 Tesla, ©o is in MHz range, which is the radio frequency band of the 
electromagnetic spectrum. The brief application of Hi rotating at ®o is called RF pulse. 

To understand the effect of the RP field Hi on the net magnetization vector M, it is 
convenient to introduce the concept of rotating frame of reference (see, Hinshaw and 
Lent [29]). This is simply a frame of reference which rotates around the static magnetic 
field Ho at the Larmour frequency coo. In the rotating frame of reference, M does not 

° . o 

appear to precess since the observer also rotates. To achieve a 90 rotation, i.e. 90 
excitation pulse, the RF amplitude Hi and duration x can be adjusted to give yHix = rc/2. 
At the end of such excitation pulse, M would be in the transverse plane and would 
appear stationary in the rotating frame of reference. In the laboratory frame of reference 
the rotation of M away from z follows a spiral trajectory over the upper half of a sphere. 
At the end of RF pulse, M will precess freely about z at its Larmor frequency. The signal 
generated by the magnetization after a RF pulse is called a Free Induction Decay (FID). 
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The free precession around z continues as the system returns to its equilibrium position. 
The behavior of M during this portion of the MR experiment is characterized by two 
processes, spin-lattice and spin-spin relaxation, which have exponential time constants Ti 
and T 2 and describe the decay and recovery of the longitudinal (z) and transverse (x and 
y) components of M, respectively. The transverse component of M precessing in an ever- 
decreasing radius circle in the transverse plane (spin-spin relaxation) while z component 
of M grows with time (spin-lattice relaxation). Because of the nature of the underlying 
process T 2 < Tj (Peters and Williams [58]). In biological tissues at typical imaging field 
strengths (0.1 - 2.0 Tesla), Ti values are of the order of 0.1-50s and T 2 is 50-1000 ms. 

It is an experimental fact that the motion of M(t) for t > to is a damped precision around 
z-axis, with the Larmor frequency, and M(t) — > Mok as t — » 00 , due to spin-lattice and 
spin-spin interactions. In 1946, F.Bloch [6] formulated the differential equations 
describing the motion of the components of the magnetization M. For the class of 
experiments in which the magnetic field is a small perturbation of a static and 
homogeneous field, i.e H(x,t) = H(t) = Hok + Hi(t), with Ho » max t |Hj(t)|, and in the 
case of M(x,t) = M(t) the following equations are basic to MR experiments and are called 
Bloch equations. 


d M x i + M j 

— M(x, t) = y M(x, t) 0 H(x, t) 

dt 1, 


(M z - M 0 )k 
T, 


(1.5.2) 


The notation, 0 in (1.5.2) refers to cross product. Componentwise, 


dt 


M x (t) = y [M y (t)H z (t) - M z (t)H y (t) - 


M x (t) 


_d 

dt 


M,(t) = y [M.COH.O) - M„(t)H,(t)] 


M,(t) 


dt 


M 


: (t) = y[M x (t)H y (t)-M y (t)H x (t)] 


M z (t) - M 0 


(1.5.3) 


Bloch equation gives a description of the time dependence of the nuclear magnetization 
M(t) in the presence of an applied magnetic field H(t). The local inhomogenities in the 
magnetic field Hok will also cause M xy (t) to decrease to zero. Denoting T 2 by T 2 * which 
also takes into account the magnetic field inhomogenities due to which the fall to zero of 
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M xy (t) will be faster, that is T 2 < T 2 (Sebastiani and Barone [66]). To incorporate this 

fact in the Bloch equations, let W(x,t) be the vector with components 

"_ M x (x,t) M y (x,t) M z (x,t)-M 0 (x) > 
l T 2 *(x) ’ T* (x) ’ T,(x) )■ 

We assume that relaxation’s due to spin-lattice and spin-spin interactions and magnetic 
field imperfections can be taken into account by adding to the right-hand side of (1.5.2) 
the term W(x,t). Then the Block equations become, 
d. 

— M(x,t) = y M(x,t) ® H(x,t) + W(x,t) (1.5.4) 

where 

H(x,t) = Hok + Hi(x, t) with H 0 » max Xjt | Hi(x,t)|. These Bloch equations provide us 
with a macroscopic model for the net nuclear magnetization M. 

A technique known as spin echoing (see Partain et al., [56]), may be used to eliminate 
the dephasing effects, caused due to spatial inhomogenities. Here a 180-degree pulse 
applied for sometime x after the 90-degree pulse and it can be shown to reestablish the 
phase coherence x msec later (Stark et al., [69]). Immediately following the 90 degree 
RF pulse, at x = 0, magnetization is coherent (in phase) in the transverse xy plane. At 
time t = x following the 90 degree pulse, the partly dephased transverse magnetization is 
turned into mirror image positions by 180 degree RF pulse. The fast moving spins catch 
up with the slow moving spins and refocusing takes place at time t = 2x. The time 2x is 
referred to as the echo, time TE. After several repetitions of 90- 180-spin echo, the 
magnetization strength reaches steady state. M is detected at some point during the 
relaxation process of each repetition. The repetition time is usually referred to as TR. 

In 1972 Lauterbur [45] first showed that by superimposing linear field gradients on the 
main field, projections of an object could be generated from which an image can be 
reconstructed. Application of gradient alters the magnetic field in a selected direction. 
SH 

Let G = be the gradient applied in x direction, the resonance frequency becomes 

x dx 

dependent on the location of the volume element of interest with respect to the x 
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direction. Both the gradient field G and the RF field Hi are under the control of the 
experimenter. The gradient field is controlled by three input signals, called G x (t), G y (t) 
and G z (t), while the RF field is controlled by Hi(t). These signals are grouped into a four 
component vector P(t)- the pulse sequence (Partain et al., [56]). The coils which generate 
these gradients in the machine are known as the gradient coils. The ideal x-gradient coil 
causes the z-component of the magnetic field to vary linearly with x as follows, 

H(t, x) = H 0 k + G x (t) x k. 

The resonance frequency g>o is a function of position, 

©o(r, x) = y(H 0 + r G r (t)). 

The real x-gradient coil produces a field which has components in the x and y directions, 
but the uniform field in the z-direction is so strong that these components may be 
neglected. We can assume that G x (t)x is the x-gradient field produced by G x (t) the. x- 
gradient. This is true for y and z gradients also. When all the three gradient coils are 
turned on, their superimposed fields yield 

H(t, x) = H 0 k + G x (t) x k + G y (t) y k + G z (t) z k. 

The three gradients may be formally grouped into a gradient vector G(t) with components 
G x (t), G y (t) and G z (t), and we can write 

H(t, x) = (H 0 + G(t).x) k. 

The RF field Hi(t) in the rotating frame of reference is given by, 

Hi(t) = Hi [cos cot i + sin cot j], 
where Hi is the magnitude of the rotating field. 

t 

Let R(x) = | G(x, x)dx , where G(x, x) is the gradient applied in the x direction. 

0 

The net magnetic field in rotating frame of reference is 

H(x, t) = (H 0 + R) k + Hi(t) = (H 0 + R) k + Hi cos co r t i + Hi sin © r t j, 
where co r is the angular frequency of the rotating frame of reference. Substituting the net 
H which includes the effect of gradient fields, RF field in the rotating frame of reference 
and the main field Ho, in (1 .5.4), evaluating the cross product and equating the respective 
vector components we get 
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dM x 

dx 

dM y 

dx 

dM 2 

dx 


= y[M y (H 0 +R)-M Z H, sinco r t]-^h 


M, 


= y[M z H, cosco r t - M x (H 0 + R)]-— 7 


y[m x H, sinco r t-M H, cosco r t]- -^ z 


T, 


(1.5.5) 


M(t) obtained by solving the above equations under the experimental conditions will give 
the signal. The signal obtained depends on the applied pulse sequence. Consider Spin 

TE TE 

Echo (SE) sequence, 90* 180* measure, where TE is the echo time. The 

2 2 

repetition time TR, is the time taken per pulse sequence. Let Ml and Mx denote the 
longitudinal and the transverse components of M. Initially M is aligned with Ho, the 
static magnetic field, so that Ml = M and Mj = 0. Following the 90* pulse, Ml=0 and 

TE 

Mx=M. The relaxation Ti and T2 continues for a time — . 

2 

The signal intensity Ise corresponding to the SE sequence is given by (Partain et al., 
[56]), 


I SE =kM 


TR 


-(TR-TE/2)^ 


1 + e Tl -2e 


TE 


(1.5.6) 


In the above equation k is proportionality constant. The amount of signal coming from 
the volume of the sample is a function of both tissue parameters Ti and T2, and the 
experimental parameters TR and TE. M is proportional to the number of affected spins 
in each pixel, it is often called spin density (or proton density), denoted by p(x,y). 
Observed contrast in the image is due to a weighted sum of p, Tl and T2. 


Constrast = W, 


M 


AM 

J 


+ W r 


r AT, ^ 
V My 


+ W T 


A 

M 


, where W is the weighting factor. 


For a spin echo image, 
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W M = I SE = kMj 


TR \ 


1-e 


V 


J 


W T = -kM 


— — • e T ' e Tj 

vT, J 


(1.5.7) 


f 


W T = kM 


1 — e T ‘ 

v j 


tr Y te a te 


vT 2 y 


Wm is large if TR is long and TE is small. W Ti is large when TR=Ti and TE is small. 
And, W T2 is large when TR is large and TE=T 2 . By choosing a long TE and a long TR, 


the contrast in the resulting image will depend mainly on T 2 . Images constructed from 
signals, measured at a long TE and long TR are usually referred to as T 2 weighted 
images, while a short TE and short TR results in Ti weighted images, (see page 34 of this 
thesis, Image 1.5.5 Test image data set 3). 


The signal received is the sum of all individual signals from the distribution. Further, the 
observed transverse magnetization is directly proportional to the spin density, represented 
by p(x,y). MR signals obtained correspond to the whole volume of the object. Images 
corresponding to certain slices are useful in diagnosis. This amounts to isolate a slice 
parallel to the xy plane at a height zq. The slice selection is achieved by selective 
excitation (Sebastiani et al., [66]). By the application of selective pulse, the required 
cross-section can be selectively excited. The received signal corresponds to the selected 
slice. In MR imaging the main objective is to discriminate the variations of spin density 
within the selected slice as a function of position in the transverse plane. Frequency 
encoding is used to determine the variation of p(x,y) along x axis. Turning the x gradient 
after excitation, during the spin echo does frequency-encoding (Stark et al., [69]). 
Application of x gradient breaks the selected slice into lines of spins transverse to the x- 
axis. Each line is comprised of all spins whose x coordinates have a given value. Images 
can be obtained, from frequency encoded data by using projection reconstruction 
methods developed by Lauterbuer [45]. However, a most of the modem imaging systems 
uses a technique called spin warp imaging (see Mansfield et al., [50]). Here phase 
encoding is used to determine the variation of p(x,y) along y axis. Phase encoding of the 
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y-axis requires many echoes to be acquired. The y gradient is turned on for some time t 
at a strength, which varies from one echo to another. During the phase encoding stage, 
the lines of spins corresponding to different values of the y coordinate are spinning at 
different rates, due to the applied y gradient. At the time when the y gradient is turned 
off, the signal received corresponds to varying frequencies in x direction and varying 
phases in y-direction of the spins. 

Consider 2D distribution of spins, the total signal S(t), is the sum of all individual signals 
from the distribution, as represented by, 

s(t) = JJl SE (x,y)e T 2 - e -i(x.k x (o +y .k y( t)) dxdy; (1.5. 8) 

xy. 

where Ise given by (1.5.6), is valid for each point of the cross-section. 

t t 

In (1.5.8), k x (t) = Y J G x (r)dr and k y (t) = y{G y (x)dx. 

0 0 

We note that the MRI raw signal is the Fourier transform of the object evaluated at the 
Fourier space (k-space) coordinates k x (t) and k y (t). Signals are being measured while 
gradients are applied, the number of phase cycles multiplying the spin distribution is 
constantly changing. Therefore, at each instant, the measured signal represents a single 
sample of the 2D Fourier transform on k-space. By appropriately designing the pulse 
sequence (see Partain et al., [56]), these samples can be made to cover all of k-space. 



Raw data and Image 


Image 1.5.1 




Image 1.5.2 


Image 1.5.1 is the magnitude of the raw data and 1.5.2 is the magnitude image obtained 
from the raw data. We see that in the raw data maximum information is around the 
center, i.e. low frequency components have higher magnitude than the high frequency 
components. Based on the frequency components, a variable length coding (VLC) could 
be devised so that the low frequency components are stored at a higher precision than the 
high frequency components, resulting in compression. 

1.6 Experimental Setup: 

All the experiments presented in this thesis are done on the actual raw data and images 
obtained from MRI machine, at Sanjay Gandhi Post Graduate Institute of Medical 
Sciences (S.G.P.G.I.M.S), Lucknow, India. For most clinical applications, 256 x 256 
images are used and our experiments are performed on 256 x 256 image data. We use 
three different image data sets for three different applications [cf: 1.8 (cl),(c2),(c3) ]. We 
use 2 byte images (depth, 16 bits per pixel) for images to be used for quantitative 
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estimations and 1 byte image (depth, 8 bits per pixel) for images to be used for visual 
inspection. The performance of certain standard compression methods is tested on the 
test image data. Three types of test images are used for experiments in the thesis. Test 
image set 1 consist of magnitude reconstructions corresponding to four different cross 
sections (image 1.5.3). Images in set 2 (image 1.5.4) are obtained following the general 
methodology proposed by Rathore [62], [37] for sharpening images to improve the visual 
diagnostic resolution. In the test image set III (image 1.5.5), we consider T), T 2 and spin 
density weighted images. These images corresponding to the same cross sections are 
considered. These images are obtained by using suitable experimental parameters TE, 
TR. The TE, TR values corresponding to each of these images are given below the 
images. Performances of standard compression routines are tested for the test images and 
are shown in tables (1.5.1)-(1.5.3). Huffman coding and Arithmetic coding are 
techniques which reduces the statistical redundancies considering the data set as one 
arising from a zero order Markov process. Therefore, the average bit rate achieved 
through Huffman and Arithmetic coding cannot be less than the entropy of the image 
data. We see from tables (1.5.1)-(1.5.3) that bit rate achieved by Huffman and 
Arithmetic coding are quite close to the entropy. JPEG uses transform-coding techniques 
based on DCT. In case of JPEG compression, redundancy measures based on 
smoothness (defined more precisely in chapter 3, of the present thesis) are more 
appropriate. For smooth images we observe the compression is more in case of JPEG. 
Another compression method, which is widely used, for compression of gray level image 
is GIF (cf: 1.4.2). It is based on dictionary based coding techniques. Higher the 
correlation between adjacent pixels, i.e. higher the smoothness of images, methods like 
JPEG and GIF attains higher compression ratios. The smoothness measures introduced 
in chapter 3 are consistent with the performance of GIF. 

1.7 Thesis organization: 

Rest of the thesis is organized as follows: 

Chapter II: is a study of compression of MRI images by a sequence of trigonometric 
polynomials arising from the Fourier series of the underlying analog intensity distribution 
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function. The first section concerns with certain preliminaries including the definition of 
Besov Spaces. These spaces are identified with a weighted sequence space l“(X) , their 

norm equivalence is obtained in section two. In section three, we prove a Hardy type 
inequality and some lemmas including bound on Fourier coefficients by integral modulus 
of smoothness. These lemmas are used in our main theorem and subsequent storage 
bounds. Direct and inverse theorems which relates MR images in Besov spaces and their 
degree of approximation by partial sums of Fourier series are presented in section 4. For 
the case 1 < p < oo, approximation by partial sums of Fourier series are used and for the 
case p = 1 , oo Steklov approximations are used. Compression algorithm and theoretical 
storage bounds for MR images are obtained in section five. Section six deals with the 
experimental implementation of the proposed compression algorithm. In view of some 
advantages by storing the raw signal data instead of images, we also study the storage of 
MRI raw data. 

Chapter III: is a study of compression methods using finite differences. For sufficiently 
smooth images, higher compression ratios are obtainable through finite difference 
methods. We store k-th differences instead of pixel values, the value of k is determined 
by the smoothness of the image. In the first section, the method of storing k-th 
differences instead of pixel values for compression is discussed. For images belonging to 
the space X k,a we estimate the smoothness by measuring the parameter a. Estimated 
smoothness of test images, compression algorithm and the corresponding storage 
requirements are presented in section two. In section three an inductive proof of a lemma 
due to S.N.Bemstein and certain other inequalities are presented. The reconstruction 
procedure adopted is approximation by local interpolation by k-th degree algebraic 
polynomials. Direct and inverse theorems for the images arising from MRI are obtained 
in section four. In section five, the experimental implementation of the proposed 
compression algorithm is carried out on the test images. 

Chapter IV: is a study of compression by reversible smoothing. Four different 
approximation process Fejer, Poison, Jackson and Steklov means are used for smoothing 
the original images. The window functions and the multipliers corresponding to the four 
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approximation processes are presented in the first section. In section two, compression 
algorithm and the performance of the approximation processes at different level of 
smoothing are discussed; they are followed by experimental implementations of the 
proposed compression strategy. In section three, we obtain a relation connecting depth, 
resolution and smoothness of a digital image. In section four we obtain storage bound for 
the compression of images through wavelets. Finally we experimentally compare the 
performance of wavelet compression methods and the compression strategy based on 
different approximation processes. 



Magnitude reconstructions -Uncompressed (Test images set 1) 




68.img uncompressed 


72.img uncompressed 


64.img uncompressed 


60.img uncompressed 


Image 1.5.3 






Sharpened images (Test images set 2) 




img.36 


img.3612 


img.36 123 


img.36 1234 


Image 1.5.4 






PD T1 and T2 Weighted Images (Test image set 3) 



PD (TE=12 msec, TR=2200 msec) 


T1 (TE=14 msec, TR=1000 msec) 



T2 (TE=80 msec, TR=2200 msec) 


Image 1.5.5 







Table 1.5.1 




TEST IMAGE SET 2 



Table 1.5.3 





CHAPTER II: COMPRESSION THROUGH 
TRIGONOMETRIC POLYNOMIALS 


2.0 Introduction: 

A study of Besov Spaces B“ ,k (L p (I)) (defined more precisely in section 1 of this chapter) 

has been made in the context of compression of wavelet decompositions by R.A.DeVore, 
BJawerth and V.A.Popov [17] in 1992. They have characterized functions having 
certain prescribed degree of approximation by dyadic splines. The interpolation spaces 
between a pair of Besov Spaces have also been studied by R.A.DeVore and V.A.Popov 
[16]. In the present chapter, besides characterizing MRI images by a sequence of 
trigonometric polynomials obtained from Fourier Series of the cross-section f, we obtain 
theoretical storage bounds for the proposed compression algorithm. 

The remainder of this chapter is organized as follows. The first section concerns with 
certain preliminaries including the definition of Besov Spaces. These spaces are 
identified with a weighted sequence space 1“ (X) (defined more precisely in section 2), 

their norm equivalence is obtained in section two. In section three, we prove a Hardy 
type inequality and some lemmas including bound on Fourier coefficients by integral 
modulus of smoothness. These lemmas are used in our main theorem and subsequent 
storage bounds. Direct and inverse theorems, which relate MR images in Besov spaces 
and their degree of approximation by partial sums of Fourier series, are presented in 
section four. Compression algorithm and theoretical storage bounds for MR images are 
obtained in section five. Section six deals with the compression of the raw data and the 
chapter ends with section seven where the experiments and results are discussed. 
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We use the following notations throughout this chapter. The domain under consideration 
is the square, I=[0,27t] 2 . Let f be the analog intensity distribution function of the 
underlying cross section. It is viewed as a mapping from I to the interval [0,1], 

f : [0,2 tc] 2 -> [0,1]. 

The function f(x) is of two variables, with x = (xi,x 2 ) g I. Let n = (ni,n 2 ) g Z\ and h = 
(hi, h 2 ) g 91 2 . The operations addition and dot product are given by, x + h = (xi+hi, 
x 2 +h 2 ), x.h = xih) + x 2 h 2 , x.n = X[ni+x 2 n 2 , for t >0, |h| < t => |hj| < t, |h 2 | < t and for any 
positive integer k, kh = (khi , kh 2 ). 


2.1 Besov Spaces: 

For a fixed a > 0 and an integer k > 0. The k-th forward difference with a parameter 
h=(hi,h 2 ) € I is given by 

A° h f(x):=f(x), 

andforj = 1, ,k, A J h f(x):= A ]]" 1 f(x + h)-A j h _1 fix), where x = (x'i,x 2 ) e I. 

For 1 < p < oo, the k-th modulus of smoothness in Lp(I) is defined as 


and for p = oo, 


co k( f ’ t ) P =sup f| A k f(x)] p dx 

(hist 1 f . 


© k (f,t)=sup(|A k h f(x)|) 

ih|<t 


1 


( 2 . 1 . 1 ) 


( 2 . 1 . 2 ) 


The Besov space B“' k (L (I)) (see, R.A.DeVore, etal.,[17]) norm is defined as 


B.‘ k (L p (I)) 


' 2 k jj 

H|f|lp+ -J[t— €D k Cf,t) p f — 

L n ^ 


for 1 < p < oo, 1 < q < oo and for q = oo, 

i! f n B «^(L (i» : =ll f Hp +sup[t" a co k (f,t)p] 


t>0 


f G B“’ k (L p (I)) if it f !I B “- k (L p (i)) <a) 


(2.1.3) 


(2.1.4) 
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The limits of integral, in the right hand side of (2.1.3) go from 0 to co in the definition 
used by R.A. DeVore [17]. Since we are dealing with 2 ti periodic functions we let the 
limits to go from 0 to 2n. Although Besov spaces are also defined for 0 < q < 1, we 
restrict ourselves to the normed spaces, i.e. q > 1. We note that when q = oo, we require 
<a k (f,t) p to decay at least as fast as 0(t“) as t— >0, whereas when q < oo we require 
co k (f, t) p to decay at a slightly faster rate. 

The spaces B“ ,k (L p (I)) and B“ ,k '(L p (I)) are related as follows, 

® k+ i(f,t) p <2co k (f,t) p 

fe B“’ k (L p (I)) =>fe B“’ k+1 (L p (I)). (2.1.5) 

It follows by induction that for any k' > k, f will be in B“ ,k ' (L p (I)) if f is in B“ ,k (L p (I)) . If 
both k and k' are strictly greater than a, then B“ ,k (L p (I)) is equivalent to B“ ,k ’ (L p (I)) . 
This result is presented in the form of a theorem. 


Theorem 2.1.1: If k' > k > a then 


M| Ee k (L p (D) 


To prove this theorem we need three lemmas, which are as follows. 


Lemma 2.1.1: The following identity for k-th differences of a function of two variable f, 
holds 

k 


A‘ b f<x)=y 




A k f(x + oh). 


Proof: For the case of one variable function the above identity is available (see Timan 
[71]). Here we prove the above identity for two variables by using induction on k. 

For k =1, 
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A 2h f ( x ) = f ( x i +2h,,x 2 +2h 2 )-f(x„x 2 ) 

= f(x, + 2h, ,x 2 + 2h,)-f(x, +h,,x 2 +h 2 ) + f(Xj +h 15 x 2 +h 2 )-f(x,,x 2 ) 
= A k h f(x + h) + A k f(x) 

The identity is true for k— 1 , let us assume that it is true for k-1 and we prove for k. 
A k 2h f(x) = A k ; h 'f(x + 2h)-A k 2 ;'f(x) 

= S| k 1 |[A k h ' l f( x + 2h + oh)-A k ' 1 f(x + oh)] 
o=0 V u J 


A k ’f(x + 2h + oh) - A k h ‘f(x + h + oh) + A k ! f(x + h + oh) - A k h ’f(x + oh)] 


U k 1 i[A k f(x + h + oh) + A k f(x + oh)] 

u=0 \ U J 

Krtc-iv . i ... 


= Z| K 1 ,|kf(x + uh)]+A'if(x+kh) + A , if(x)+W k ’kffx + uh)] 

o=l U V 0=1 V w ) 

nAV , A c _i\A c _]\A c N 

= 5j R f ( x + oh) + A k f (x + kh) + A k h f (x), because + 

fbT oj o-lj l u j o, 


k /M r 1 

= 2 [A k f(x + uh)J 


Lemma 2.1.2: Let 1 < p < oo, and f e Lp(I) then, 


© k (f,t) p < M k t k j|| f || p + j fl>k ^- U - p - dul, 


where the constant Mk does not depend on f. 


Proof: For the case of function of one variable and p=co, the proof is available (Timan, 
[71, pp.108]), for functions of two variables and for 1 < p < oo, the proof is similar, we 


just give the outline. 



[By Lemma 2.1.1] 
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k2 k "'co k+1 (f,h) p . 


Consider, s = (s,,e 2 ) , putting here, successively h=2 m s (m=0,l, ,k-l) and taking 

0<Si,S2 < t. We obtain system of inequalities of the form 


lA k m< , e f(x)-2 k A k 2 


f(x) < k2 k_i co k+ , (f,2 m s) . 

1 E L [0,2it-2kh] k+1V 'P 


Multiply both sides by 2’ < - rn+1 ' )k and summing for m=0,l,...,r-l and using the fact that 
C0k+i(f,u) p is non-decreasing with u, we get 


|2- rk A k r£ f(x)-A k f(x)| < k 2 t k J 


2 >k +1 (f>u) p 


which implies that 




( 2 . 1 . 6 ). 


Proceeding as in above with f(x) replaced by <|)(x) = f(27t - x), i.e. 
KXjjXj) = f(2Tc-x,,27c-x 2 ) , we obtain 


IK+wll <2- <T -‘> t ifS r +k ! t k j 


2 ',;co k+I (f,u)p 


(2.1.7) 



A k f(x)| P dx I < j|A k f(x)| P dx A k f(x)| P dx 



= Ji + J2say. 


We have. 


( 2 . 1 . 8 ) 


A k f(x)| P dx '= |A k f(27i-y)| P dy 


A k e (j)(y)| P dy . 


Hence from (2.1.6), (2.1.7), (2.1.8) and (2.1.9) it follows that 


A l 

Agf(x) < 22~ (r ~ 1)k ||f|| +2.k 2 t k f 

D ‘ J 


2 . k 2 'r®t + i( f .“)p 
* — — 


(2.1.9) 


(2.1.10) 
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*•> r\ 

We choose integer r=r(t) such that -A<2 r tk < this gives by (2.1.10) that 


Agf(x) 


L [0,2-n-ks] 


< M k t K <! 


+ 


/ K 

f 


®k+i( f 5 u ), 


u 


k+l 


-du 


( 2 . 1 . 11 ) 


Since s < t is arbitrary the result follows from (2.1.11). 

Now by using the following lemma, Hardy’s inequality (Butzer and Berens [10, pp.199]), 
in conjunction with the lemma 2. 1 .2 we will complete a proof of theorem 2.1.1. 

Lemma 2,1.3 (Hardy’s inequality'): Let a > 0, 1 < q < oo. If v|/(s) is a non-negative 
measurable function on (0,°o) (measurable with respect to the measure ds/s), then 


and 


o v. o 


00/ 00 


]R<) f' s 


q ds ( 
s ( 


a r / x ds ] dt 

'Kt 


0 V t 




(2.1.12) 


(2.1.13) 


holds. 


Proof of Theorem 2.1.1: 


co k+1 (f,t) p <2co k (f,t) 

I f IL.™ <C||flU, 




Conversely, by Lemma 2.1.2, 


2tc ij ] q 271 

[(t-XCf.t),) 1 - s j 

n t n 


r“ 


2n/2k 


M k t 


®k +! ( f > u ) f 


u 


k+l 


-du 


p 

1 

q 'N 


dt 



t 


J 


ail 

<M k 2 ( > 


lJ t (k-a )q -. dt+Mk | t k-a |°>k t l(p- U l ^ 


dt 

t 


\r 
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-^k,a f p +M k I, 


(2.1.14) 


where, I,= 

^2rc 

j 

^.k-a r®k+l(^>^0 

dt 


1° 

' u k u 

t 

J 


We evaluate Ii by using Lemma 2.1.3. We define 


V(s) = 


I.= 


V ° L 


J t 1 '” Jv(s) 


ds 

s 


scj J(s l ->(8)f- 


t 

ds 
s 


|s' k © k+1 (f,s) p 0 < S < 271 

0 s > 2n 


[By Lemma 2.1.3] 
[Where the constant C a ,k depends only on a and k] 


' \ q f 2n i 

C a , k ( }(s k - a s- k 0) k+1 (f,s) p ) q ^| = C ak [ }(s-“co k+1 (f,s) p ) q ^ 


(2.1.15) 


Substituting (2. 1 . 1 5) in (2. 1 . 14) we get, 

H f ilB“' kl (L p (I)) - C 1 1 1 f l ! B“* k+ ' (L p (I)) 

Hence the proof. 

In the next section we define certain sequence norm and will show its equivalence with 
this Besov space norm. 


2.2 Equivalent Norms: 

In the present section we identify Besov Space B“(L p (I)) norm with a weighted sequence 
space 1 “ norm by proving their norm equivalence. 
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Definition 2.2.1: 

If a = {a n } is a sequence whose component functions are in the quasi normed space X, we 
use the 1 “ (X) norms, for 1 < q < oo 


hu- sKwixr^r p- 2 - 1 ) 

v n =i lL J 

with the usual change to a supremum'norm when q = oo. When {a n } is a sequence of real 
numbers, we replace ||a n ||x by |a n | in (2.2.1) and denote the resulting norm by ||(a n )|| „ . 

Theorem 2.2.1: For f e Lp(I), 

f e B“’ r (L (I)) if and only if jcoff ,— 1 1 e 1 “ 


Proof: 


Let ^co f, — > e 1“ . Therefore we have V n a co r f, — — <oo. 

1 I n J pJ „ V n JpJ n 


2*/ 

oo /k . 


J[t" a ® r (f,t) p ] q — =£ I [r a co r (f,t) p ] q - 

0 1 k=l 2 k / 1 

/k+1 

1 °° ( 1 

<7—7 — V (k + l) aq © f,— - 

(2*rtr l kj p k 


[o) r non decreasing] 


•• 2 T * < , 2 ” ( k+1 ) _ 1 

’ 2} t _ (k + l)k 2% k 


00 f 1 

<C a , q X k a co r f, — — <00 by hypothesis. 

k=i j_ v k J p J k 

Therefore it follows that f e B“(L P (I)) . 


Conversely, let f € B“(L (I)) , 


2 V 

dt ^ ?| 


i t n = 1 Oir / 


^ ^ a ( - 2tt ) 1 

> 7 n o. f, — — n 

. l n + ljp n(n + 1) 
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2 aq 


■z 

n=l 


(n + l) a co r 


. 271 

l ’n + U 


1 

(n + 1) 


Therefore, 


f , 

(. 2n 

,1 


f,— 


l 

l n ) 

P J 


ei :- 


Hence the proof. 


In the next section we shall prove some lemmas which will be used along with this norm 
equivalence to characterize Magnetic Resonance Images in Besov spaces. 


2.3 Hardy type inequality and bound on Fourier Coefficients: 


In this section we shall obtain a bound for Fourier coefficients in terms of modulus of 
smoothness and in what follows, we shall prove Hardy type inequality which will be used 
in proving the inverse theorem in the next section. 

Lemma 2.3.1: Let c > 1, p >1 and S n =a, +a 2 + --- + a n , where aj>0then 

fV c S£ < K p>c Jn- C (na n ) p , 

n=l n=l 

where K p>c is a constant depending only on p and c. 

Proof : 

Let So = 0. Define <j) n = n' c + (n+l)' c + then <t> n < k n’" c . This follows from the fact 

that n" c +(n + l) _c + < f— = — — we take k = — — . 

y 1 J X c 1 - C 1 - c 

n 

m m 

£n-‘S p = ^(♦ n -* n+I )sS 
1 1 

m m 

1 1 

m m+1 

12 * 

m 

* ♦,sr+S+.(s:-s:.,) 

2 
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m 

£ *,(sr-s;)+l>.(s:-s:_ 1 ) 

2 

m 

1 

m 

1 

111 

1 

m 

- k P,cE n '" Cs r la .. where k p>c = -£- 

1 l ~ C 

m _£ JL 

- k P,cS n Pn p s n"'( na n)> where - + — = 1 

i P P 


= k p ,£n P K)n P 'S;' 




( ^ 

n 

^ m 'N 

- k p,C 

E n ’ C ( na n) P 

F 

In-s; 


v 1 ) 


\ 1 ) 



( m 

> 

D 

f m 

A 

=> 

E n_C S P j 

^ k p,c 

E n-‘(naj 


V I 

) 


V i 

J 


m m / \ p 

E n ' CS ^ K P ,cE n ' 1 na ») P > where K p , c =(-£-] 

1 1 VIC/ 


Hence the proof. 


[By Holders inequality] 


Corollary 2.3.1: 

n n 

Let c > 1, q >1 and S n = ^^a ;j _ where ay > 0 then 

i=i j=i ■ 

IX C S n ~^q, C X n C ( n 2 a nn ) q » 

n=l n=l 

where K p c is a constant depending only on q and c. 
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Proof: 

Let S n = P" + • • • + P n n , where P" = J a kj 

k=l 

Jn" c S^ ^K qc ^n' c (nP n n ) q , applying lemma 2.3.1 

n=l n=l 

= K q.c £ > where T n = J i 

n=l 

<K^X„-(n J a m r 


na 


k=l 


n=l 


The following lemma bounds the Fourier coefficients. 

Lemma 2.3.2: For f e Lp(I) , 1 < p < oo , and a given k e Z + we have for n =£ 0 

|f(n)| < C k co k (f,t) , 


where n = (ni, a 2 ), t = 


independent of f. 


n 


2(1 n, | + | n 2 |) 


> 0. The constant Ck depends only on k and is 


Proof: 

Let us consider, I = [0,2 ti] 2 , h = (hi, I 12 ) e I, n = (nj, n 2 ) e , x = (xj, X 2 ) € I, 
dx=dxidx 2 , and the dot product n.x = niXi + n 2 X 2 . We have, 

(a* f)'(n)=-L f(A k „-'f(x + h) - A k h "'f(x))e"“ x dx 

471 j * 

= J_ e inh [A k h -'f(x)e~ inx dx l — fA k -'f(x)e- inx dx 

Arr 2 A'TT 1 


=J— (e in h -l)fA k -'f(x)e- in x dx 
4i t , J 

= _L( e in h -i) k Jf(x)e" in x dx 


•h f in.h in. 


e 2 — e 2 


■i 2?t 

«") J Jf|A‘ f < x )|dx 

10 
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2 sin 


n.h 


f( n )|^— r ( 4 ^ 2 >'[ J|A k h f(x)|' 


[By, Holder’s inequality] 


Let t = 


n 


i,i -r, and h = (h„h 2 ) = 
2(| n, | + | n 2 |) 

take hi = -hi and if n 2 < 0 , h 2 = -h 2 . 


n 


71 


2(n, +n 2 ) 2(n, +n 2 ) y 


If ni < 0, then we 


2 k 


f(n) < G p co k (f,t), where C p = 


(471 2 ) 


f( n )| ^ Cp.k^klfn), where C p>k = — r 

2 k (4:r 2 )r 


Hence the proof. 

In the next section we shall present the direct and inverse theorem. 


2.4 Direct and Inverse Theorems: 


Here for the direct theorem we deal the cases 1 < p < 00 and p = 1, 00 separately. For the 
case 1< p < co we use partial sums of Fourier series and for p =1, 00 Steklov 
approximations are used to establish the results. Inverse results are obtained for the best 
approximation trigonometric polynomials. To prove our main theorem, certain results 
about trigonometric polynomials of best approximation for functions of two variables are 
required, and these are available in the literature. Here we just state the results. These 
results are based on partial modulus of smoothness (Timan, [71]). For functions of two 
variables the partial modulus of smoothness is defined as follows. 


Definition 2.4.1 (Partial Modulus of Smoothness): 


<n k (f;u,0) = sup sup 

E(-D k -" 

f k ] 

f(x + ut,y) 

y 

11=0 



-- 

k 

f k l 


6 k (f ;0, v) = sup sup 

Eh) 1 "” 

f(x,y + ot) 

x |t|£v 

0=0 

l u J 



and for 1 < p < 00, 
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k lk -/^ 


®k(f;u,o)p =su P | 


®k( f ;°> v )p =su P j52(-i) k_u 

l>A i u=0 


fk\ 


v u y 


f(x, + ut,x 2 ) 


f(x, +ot,x 2 ) 


dx 


dx 


Definition 2.4.2 (Total Modulus of Smoothness): 

f(x, +uh,,x 2 +uh 2 ) 


© k (f;u, v) = sup sup 

|hj|<u |h 2 |<v| 

and for 1 < p < oo ? 


k fh-\ 

E(-D k - J 

u=0 


v u y 


J jA 

© k (f;u, v )p = sup sup J^](-l) k ' u f(x, +uh 1 ,x 2 +uh 2 )dx 

|h,|<u|h 2 |<v ' u=0 yOy 


Definition 2.4.3 (Mixed Modulus of Smoothness): 


®k,i( f ; u > v ) = sup sup A h A h f(x,,x 2 ) 

|hj|<u |h 2 |<v * 1 


= sup sup 

|h t |<u |h 2 (<v| 


k 1 


u=0 ji=0 


m 


f (x^ + oh 1 ,x 2 + oh 2 ) 


and for 1 < p < oo, 


k 1 

|h, |<u |h 2 |Sv |" ' o=0 M=0 


®k,i( f ; u A)p = sup sup j i2(-i) k+1 " u ' 11 


rkYi 






f(Xj +oh,,x 2 +uh 2 )dx 


The relation between moduli of smoothness is as follows (Timan, [71]). 

k f’k') 

max{© k (f;u,0), © k (f;0,v)} < © k (f;u,v) < ^ © k . u , u (f;u,v) (2.4.1) 

o=0 \PJ 

For the case when u = v = t, we have cq.(f;u,v)=cq.(f;t), where © k (f ; t) for functions of two 
variable is defined by equations (2.1.1) and (2.1.2). Let E*(f) p =||f-T n || p where the 

function f is of two variables, and T n (x) is the trigonometric polynomial of best 
approximation, or order n,n in the corresponding variables xi,x 2 . 


f& t x W in-nfbTTsr WT 


V o«g JSajlfeft 

VtfW ASM 
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Theorem 2.4.1: Consider the function f(x) with x = (xj, x 2 ), defined throughout the 
whole space of the variables xi and x 2 and having period 2 n in each of these variables. If 
1 < p < oo, and if f e Lp(I), then for any integers kj, k 2 , n > 0, with k = (ki, k 2 ) we have 


El(f) p ^C k 


on 


v 


,o 

n + l 


+ ©k. 


1 ^ 

f;0,-4 

n + lj 


A 


'P V “ ■ J 

where Ck is a constant depending only on k, and ct)k(f; ti,t 2 ) p is the corresponding integral 
partial modulus of smoothness. 


From the above theorem and (2.4.1) it is clear that, for 1 < p < oo, 

. E* n (f) p <C k (o k (f,t) p , (2.4.2) 

where t = — - — and k = ki = k 2 . The <»k(f ; t) in the right hand side of (2.4.2) is the 
n + l 

modulus of smoothness for functions of two variables, as defined in (2.1.1). Next we 
state a theorem which relates approximation by partial sums of Fourier series and the best 
approximation in case of 1 < q < oo. 

Theorem 2.4.2: If 1 < q < oo and if the function f(x), x = (xi, x 2 ), periodic in both the 
variables of period 2u belongs to L q on I, then for any n > 0, 

l 

| J|f(x> - s; (f;x)|'dx I" < c, e ; (f>, , 

where C q is a constant depending only on q and S* (f ; x) is the partial sum of order n,n of 
the Fourier Series of f(xi,x 2 ). 

Proof: 

Let T n (x) be the trigonometric polynomial of best approximation, of order n,n in the 
corresponding variables xi,x 2 . We have, 



Chapter 2 


53 


IM;mlL =l f -T. +T » -s;,(f)|| q <(i + |s;|J|f-T,| 

s C 1 ||f-T„|l S* is a bounded operator in L q (I) for 1 < q < °o 

and the constant C q depends only on q 

- C q E‘(f) q 


Next we state the converse theorem for best approximations by trigonometric 
polynomials. 


Theorem 2.4.3: If for a certain q (1 < q < oo) the function f(x) belongs to the class L q 
over the square I, and E* (f) (n=0,l,2 ) is the sequence of best approximation by 

trigonometric polynomials T n n2 (x), of order ni,n 2 in the corresponding variables xi and 
X 2 , then whatever the integers ki,k 2 maybe, the inequality 


co 


k, t k 2 


1 1 


P n n 

XT' V 


CAN 

V. .13 31 J n n o, =0u,=0 


u 1= 0u 2 =0 


f in 

remains valid, where C k t is a constant depending only on ki, k 2 and cb k k f,— 

1,1 " V n nj 


is the corresponding mixed modulus of smoothness. 


Theorem 2.4.4: If f € B“' k (L p (I)) for a certain p (1 < p < oo), q (1 < q < co) and for any k 

and a such that 0 < a < k, then |E n (f ) p } n el“ where En(f) p = || f - S n (f) || p and S n (f) is the 

partial sum of order n, n in the corresponding variables x l5 X 2 , of the Fourier series for the 
function f. 

Proof: 


£|n*E.<0,N S C,XkE; ( f)Ji 

n=l n n=l 11 


OO 


"f,- 



n co k 


n=l 


k n y 



n 


[By theorem 2.4.2] 


[By (2.4.2)] 


< OO 


[By theorem 2.2.1] 
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=> |E n (f ) p } n € 1“ . Hence the proof. 

For the case p = l,oo we make use of Steklov approximations 


Definition 2.4.4 (W m p (I) Space): 

Let 1 < p < co, W ni,p (I) is space of all functions for which ||f||m, p < °o, where the norm is 
given by, 

ML 

Steklov Approximations: 

Definition 2.4.5: 

The m-th Steklov approximations {f n , m }, (r) > 0, m = 1,2. . .) is defined as 

— — f A 

^ T) > \ -m m 111 

- J. • J fW-(-ir a" f(x) d n , •dn„. 

v J I 

Fori < p < oo, for every f e Lp(I), {f^m }eW m,p (I) and we have the following: 
Theorem 2.4.5: For every f € Lp(I) 1 < p < co and rj > 0, 

|| f — f n.m| p - CO m(f> r l)p- 

Proof : 


fn.m( X )-f(x) = 



*1 

JL 


f. 

m 

f 

Vm, 

) J 

• 0 

J 

0 


Using Generalized Minkowski’s inequality we have, 


f_„-f <1- 


nYlnr 


m; ymj 


* sup ||A:f(x)| 


A" f(x) 
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Theorem 2.4.6: If f e B“’ k (L p (I)) for a certain p (1 < p < °o), q (1 < q < co) and for any k 
and a such that 0 < a < k, then |e„(f) p } n e 1“ where e n (f) P = || f - f n || p and f n ’s 

k 

correspond to f n ,k, the k-th Steklov approximation of the function f, for q = — . 

n 

Proof: 

The proof of theorem 2.4.6 is similar to that of theorem 2.4.4. The proof for this case 
uses theorem 2.4.5. 


Now we shall state and prove the inverse theorem. Inverse theorems are proved for the 
best approximation trigonometric polynomials. 

Theorem 2.4.7: If for a certain p (1 < p < co) the function f e Lp(I) and {e* (f) p } n el“ 

where E p (f) p =.||f-T n || and T n (x) with x=(x u x 2 ) is the best approximation 

trigonometric polynomial, of order n,n in the corresponding variables xi,x 2 , for the 
function f, then for all k, and a such that 0 < a < k, the function f eB“’ k (L p (I)) . 


Proof : 

00 

If , we show that, ^ 


n=l 



f,-' 


Ot - 

n ® k 



l n; 

P _ 


— < oo , then the above theorem follows from 
n 


the theorem 2.2.1. 

Consider, eo k and <» ku the total and mixed modulus of smoothness respectively (cf: 
definitions 2.4.2 and 2.4.3). For n > 0, we have, 


f 


<o v 


1 


f;- =o k 
V n ) 


1 1 


^fk 


V n nj 


G) 


k-u,u 


f .II 

i 5 > 

V n n) 


‘f k 


- ZL |Hr Z Z (u > + 1)k ~ u-1 + 1)0-1 (f) £ By theorem 2A3) 

n. o 1 =o\) 2 -o 


u=o 




u=0 




n n 




ft Ui=0u,=0 
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i2 ‘ k rfSE" l ' 2 E' (0 


u l =0 u 2 =0 


[where E UiUj (f) is the best approximation of order 1 ) 1,02 in the 


Therefore, we have 


corresponding variables xi,X 2 of f(x).] 


Z 

n=l 


n co, 



V n 



Z 

n-1 


C 


"7EZ"‘X,(f), 

11 n. =1 it. =1 


Uj — 1 u 2 =l 


n 


< Ck^n^-'SJ 

n=l 


[whereS 0 = ]T Ja ;j and a ;j = n k “ 2 E*.(f) p ] 

i=l j=l 

oo 

< C k J]n _(I+(k " a)<!) S2 

n=l 

[(1 + (k-a)q ) > 1 v q > 0 and 0 < a < k] 

Applying the corollary of lemma 2.3 . 1 we get, 

-Cq,kZ n " (1+(ka>q> ( n2a ™) q 

n=l 

< c ,, k y[n“E;(f) t ] q n-'n- k ’n k ' 

n=l 


^zk E ; ( f) t r 


n-1 


l 

n 


< oo, by hypothesis. 

Here the constant C q> k depends only on q and k and is independent of the function f. 
=> f e B“' k (L p (I)) and hence the proof. 


In the next section we propose a compression algorithm for MR images and we obtain 
storage bounds. 
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2.5 Compression Algorithm and Storage bounds: 

Let the given MRI image be 

f(x,y) - ( 2 . 5 . 1 ) 

We apply an invertible transformation on the Fourier coefficient’s, to get a new function 
f . This transformation could be based on a known convolving function or otherwise, so 
that 

f(x,y) = 2> i *e' u ” ky> < 2 - 5 - 2 ) 

j.k 

with C jjk | < |c j k | holding. The next step in the compression algorithm is to modify c j k ’s 

to c* k ’s such that only non zero c* k ’s are to be retained. Normally the zero coefficients 

occur in clusters so that the storage indexing is a minimal.The Compression Ratio (C.R) 
is given by 

Uncompressed image size-Compressed image size 5 3) 

Uncompressed image size 

A higher compression ratio in the method is obtained indeed because of the number of 
vanishing c* k ’s . The compressed image data file contains a header. Patient information 

and certain parameters are stored in the header. We store the zero frequency term in the 
header and all other frequency terms are stored by using the following algorithm. 

Algorithm 2.5.1 : 

Step 1 : Calculate c Jik ’s from the pixel data p Jik ’s using FFT. 

Step 2 : Calculat c J k ’s from Cj, k ’s for the chosen convolving function. 

Step 3 : Choose a positive integer N and numbers 1 < p < 00, a > 0. Choose q such that 
ocq > 2 and quantized coefficients c* k that satisfy 
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(iii)c jk = Oif |c jik |<l 

IN 

Step 4: Store non zero c* k ’s by a variable length coding procedure (cf: 1 .2). 

Now we give a bound on storage space required for the above algorithm, here the cj k ’s 

are the Fourier coefficients of the pixel data and f is the approximation of f by partial 
sums of Fourier series. 

Definition 2.5.1: 

Let n = (ni, n 2 ) 

C(|n|) = max{f(n 1 ,n 2 ) a |n, | + |n 2 |=n}. 

From the above definition it is clear that, 

E E|f(n„n,)| £ 241 n | f, m (|n |) + f,”(0), (2.5.4) 

n x eZn 2 eZ jn|=l 

/v /s 

where f s m (0) is the maximum of the set containing one element, f (0) . It is stored in the 

header. We do not include the zero frequency information, in our compression estimates. 
A small storage space is usually reserved for the header information. 

By lemma 2.3.2, we have, 

f"(| n |)<C lffik ff;-t-l . (2.5.5) 

v I n \ J 

The above relation is used in obtaining storage bounds for the algorithm 2.5.1 


Theorem 2.5.1 (Storage bound): 

If f.e B“ ,k (L p (I)) fori <p<oo, 0<q<oo, 0<a<k and aq > 2, then r\ the number of 
non zero coefficients c* k ‘s in the above algorithm satisfies, 
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where ! f l ! ,;-, L ,, l)) =[][<'X(f,t) ( I' 

Proof: 


dt 

t 


1 

q 


Let N be the quantization parameter in the algorithm 2.5.1 and x\ is the n um ber of non- 
zero coefficients. Let the zero frequency information f(0) , be stored in the header as per 
our compression procedure. We have 


-H< 

N 


Z 




< 


c j.k 


II' 

j.k 


•j.k 


n { eZn 2 eZ 


|n|=l 


f "Cl n |) 


[■-• |c* k | < Jcj- k | by algorithm 2.5.1 

c j, k =f(j,k)] 

[By (2.5.4)] 


[ For | n |= 0, storage will be in the header] 


- C^4 1 n | co, 

|n|=l 

<C^| n | aq_1 
N-i 


V" 

V I n \J 


a>. 


v f: |»U 


*<2 

l»l=l 


n r CO, 


f 1 ' 


fr 


n \J 


n 


r\ * C k N|f| q B „., 


(L P (I» 




[By (2.5.5)] 

[•/ aq > 2 => aq-1 > 1] 


Poisson Approximation Process: 

Poisson kernel (see Zygmund [79]), is convolved with the pixel data to obtain Poisson’s 
approximation for the image. This is achieved by multiplying the Fourier coefficients 
f(m,n) by r |ra|+|n| , for 0 < r < 1. This convolution imposes smoothness in the image. 
Lower the value of r higher the smoothness. The value of r is usually taken to be near 1 
to avoid blurring. Algorithm 2.5.1 is applied for image compression with coefficients 
c j k as the Poisson coefficients, 
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** e * c m,n — f(m,n)r^ ^ . Let the Poisson filtered image be represented as, 

f(x,y) = £c. k e i(jx+ky) . 

Now the quantized coefficients c j k ‘s are chosen as per algorithm 2.5.1. 
N’ 


c j.k • 


q 1 
< 


r lj| + l k lp 
1 S.k 


| . m l , 

= c , k < — => c . k = 0 . 
I J' k l fq J- k 


The compressed image is represented as, 


f*(x,y) = Xc; k e i(jx+ky) . 

j.k 


(2.5.6) 


(2.5.7) 


Theorem 2.5.2 (Storage bound for Poisson approximation process): 

If f e B“ ,k (L p (I)) fori <p<oo, 0 < q < oo, 0 < a < k, 0<r<l, and aq > 2, then p the 

number of non zero coefficients c* k ‘s in the above algorithm satisfies, 


r| < C k N r q 


Bq ' k (L p (I)) ’ 


where 


B“' k (L„(I)) 


r ■ 


Proof: 

Let N be the quantization parameter as in the algorithm 2.5.1 and q is the number of non 
zero coefficients, we have 


JL 

N 


non zero Cj kl 




* Z. | c lkf * EM' " IlJ r,n, f(n,n 2 f ['•’ | c j.k | - | c j, k | by algorithm 2.5.4 

c j k = f(j> k) and f (0,0) is stored in header ] 


n 1 eZn 2 e 


<j4|n||r w f"(|n|) 

M=i 


V 


|n|=l 


n 


[By (2.5.3)] 


[By (2.5.4) andO <r <1] 
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SC, r’£|np-' ta k 

|ni=I 


f 1 ^ 
f; — 


< 


C, r'£ 

|n|=l 




|n| 


r 

1 \ 

M CD k 

f: 



V 

1 n ij 


[v aq > 2 => aq - 1 > 1] 


_1_ 

In 


- C k r<1 KIb“-‘(L p (!)) 

=> T1 < C k Nr q |f| q „. L(Lp(I)) . 

2.6 Compression of Raw Data: 


The usual MRI image refers to the magnitude image reconstruction from the complex 
signal data. However, MRI offers various kinds of reconstructions like .real, imaginary, 
phase, other than the usual magnitude images, and each one give information of specific 
significance for medical diagnosis. For example, the phase images (see Stark and 
Bradley [69]) help determine if a blood vessel is patent or if it is filled with clotted blood 
instead. The presence of flowing blood is suggested if a blood vessel in a phase image 
has a pattern of stripes distinct from tissues that are known to be stationary. They also 
provide information about the flow of Cerebro Spinal Fluid (CSF) (See Partain et al., 
[55]). Real images are advantageous for certain techniques notably inversion recovery 
(See Parain et al,,[56]). The information present in real, imaginary and phase images thus 
complement the strength of MR signals present in the modulus images. It is possible to 
reconstruct all these images if the raw data (MR signal) is stored instead of just the 
magnitude images. 

The raw data contains frequency information of the corresponding image. It satisfies 
decay properties, i.e. the magnitude of the data decays as the frequency increases. 
Consider concentric circles around the zero frequency term. The maximum magnitude of 
the raw data in the exterior of the circle can be shown to decay, and the corresponding 
decay could be faster if the image under investigation is sufficiently smooth. Graph 
(2.6.1), shows that for the MR images under consideration the maximum magnitude of 
the raw data in the exterior of the concentric circles centered at zero frequency, decreases 
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as the radius of the circle increases. This decay property gives rise to compression of raw 
data. In Graph 2.6. 1 Frequency circles in the x-axis refers to the radius of the concentric^ 
circles considered with center as the zero frequency term. Instead of storing the raw data 
in the matrix form, they can be stored in concentric circles and a variable length coding 
based on the apriori information obtained from graph 2.6.1 can be used so that the 
preceding zero entries corresponding to each circle need not be stored. 

Decay of Fourier Coefficients 



Frequency circles 


y - Maximum magnitude of the raw data in the exterior of frequency circle. 

Graph 2.6.1 


2.7 Experimental Results: 

Test image set 1 consists of magnitude reconstructions corresponding to four different 
cross sections. The sizes of test images are 256x256. Algorithm 2.5.1 has been 
implemented on the test image set 1. The required Fourier coefficients are computed 
using FFT routines, due to Cooley and Tuckey (see Brigham, [8]). Images obtained after 
certain compression and decompression of test images are shown in Image 2.7. 1-2. 7. 3. 
Performances of algorithm 2.5.1 are plotted in Graphs 2.7.1-2.7.4. Computation oftissue 




Chapter 2 


63 


parameters using test image set 3 are done following the method proposed by S.B.Rao et 
al., [60]. First the computations of the tissue parameters are done on the uncompressed 
images and they are compared with the parameter images computed from compressed 
images at various levels of compression. Along with the computed images (Image 2.7.4- 
2.7.7), maps of points where the computations are not carried out are given. The effect of 
compression in tissue parameter computations is listed in a table (Table 2.7.1). The test 
images used in these experiments are normalized to take two bytes per pixel. To estimate 
the amount of distortion induced because of compression decompression, we use two 
kinds of error measures, they are defined as follows. The amount of compression 
obtained is measured, by percentage compression (P.C) given by. 

Percentage Compression (P.C) = Compression Ratio (C.R) X 100 
Following quantities are used for measuring the error in decompression, 


Ei f w-f(>)i 

Relative 1 1 Error (rll ) = —— — 

all pixels 


Zifm-f©! 

& Relative 12 Error (rl2) = kl!!EEf!f ■ - 


Xi 

Eifwi 2 

all pixels J 


Eifls-fwi 

Average 11 Error (all) = . 

T otal no. of pixels 


Average 12 Error (al2) = 


El f(i)-f0)| 2 

.all pixels 


T otal no. of pixels 

In the above expressions f is the uncompressed image and f is the decompressed image. 
For the compression of Ti, T 2 and spin density weighted images, we store them in a set of 
three images corresponding to one cross section. The total size of the set is used for 
calculating the percentage compression. 



55 % Compression using Algorithm 2.5.1 



Image 2.7.1 






60 % Compression using Algorithm 2.5.1 




Image 2.7.2 






65 % Compression using Algorithm 2.5.1 



Image 2.7.3 







Percentage Compression Vs Relative Error 




Graph 2.7.2 






Graph 2.7.4 




Tl, T2 and Spin Density computed from uncompressed images 



Spin Density Map of points where computations 


are not carried out 


Image 2.7.4 







Computed from 47.38 % Compressed image set 



Map of points where computations 
are not carried out 


Image 2.7.5 







Computed from 57.68 % Compressed Image set 




Spin Density Map of points where computations 

are not carried out 


Image 2.7.6 







Computed from 66.75 % Compressed image set 



Spin Density 



Map of points where computations 
are not carried out 



Image 2.7.7 







Effect of compression in the computation of Tissue parameters 

Relative errors obtained are with respect to the parameters computed from uncompressed 
images. In the following table, P.C stands for percentage compression, rll and rl2 stands 
for relative 11 error and relative 12 error. 


P.C 

T1 

T2 1 

P.D 


rll 

rl2 

rll 

rl2 

rll 

rl2 

47.38 

0.023001 

0.108814 

0.038972 

0.082052 

0.017739 

0.110434 

57.68 

0.058826 

0.153803 

0.048279 

0.117673 

0.037774 

0.155013 

58.60 

0.070829 

0.173674 

0.048764 

0.125593 

0.040470 

0.173440 

62.41 

0.121785 

0.210169 

0.061526 

0.153636 

0.071345 

0.205007 

66.75 

0.151067 

0.259694 

0.060059 

0.158104 

0.080254 

0.237892 

71.20 

0.172699 

0.309278 

0.066849 

0.178419 

0.093853 

0.269737 


Table 2.7.1 


Compression of Raw Data 



Raw data uncompressed 


Raw data 65 % compressed 


Image 2.7.8 





Average II -error 



Raw data 70 % compressed 


Raw data 75 % compressed 


Image 2.7.9 


Percentage compression Vs Average 11 -error 

(Compression of Raw data) 



Graph 2.7.5 
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The. decompressed images at 55,60 and 65 % compression levels are shown in Image 
2.7. 1-2. 7. 3. The graphs of Percentage compression vs relative error are presented in 
Graphs 2.7.1-2.7.4. From the images and graphs we observe that there is no visual 
deterioration of images up to 60 % compression. After that the image architechture 
deteriorates slightly. Apart from test images performance of algorithm 2.5.1 has been 
tested on images corresponding to several other cross sections and the results obtained 
are similar. Upto 60 % compression is suggested for the images under consideration. 

In some specific cases of diagnosis, computed Ti,T 2 and density images are used. Three 
images (Test image set 3) which corresponds to the same cross-sections from which these 
parameters can be calulated in retrospect are compressed and stored as a set. Images 
2.7.4-2.7.7 and the corresponding error table shows that upto. 57% compression the 11- 
relative error after decompression is less than 5 %. At around 67 % compression we note 
from Table 2.7.1 that the relative 11 -error goes upto 15 %. 

Images 2.7.8-2.7.9 are reconstructed using compressed raw data. Reconstructed images 
with compressed raw data are compared with the ones obtained reconstructed from raw 
data without compression. From the following images and the error graph (Graph 2.7.5) 
we note that the distortion is not perceivable by human eye up to 65 % compression. 
Above that level, slight deterioration is observed. The test raw data permits 65 % 
compression through algorithm 2.5.1 without visual loss of information. 



CHAPTER III: COMPRESSION THROUGH 
ALGEBRAIC INTERPOLATION 


3.0 Introduction: 

The problem of image compression is viewed as a non-redundant representation of the 
underlying analog intensity distribution function, f. A digital image f , which represents 
f> is stored as samples of f called pixel values. Storage requirements of f depend on its 
resolution m, which is the size of the pixel and depth b, which is the number of distinct 
gray levels (states) the image can take. For sufficiently smooth images the number of 

states taken by the finite differences of f are less when compared to the number of states 

taken by the digital image f . Variable Length Coding (VLC) like Huffman coding [30] 
uses the peakedness of the histogram by assigning short code word for frequently 
occurring amplitudes and longer code for others. Consequently, methods based on 
VLC’s are expected to work well for the differences than on the actual images. We 
illustrate this with the help of an example. 

Consider the histogram of the image (Image 1.5.2) and the histogram of its first 
difference. Graph 3.0. 1-3. 0.2 (page-73) shows that the number of distinct states required 
for representing the first difference, is less when compared to that of the image itself. 
The method of storing first difference instead of pixel values is popularly known as 
Discrete Pulse Code Modulation (DPCM). It has been investigated since 1970 (Jain, 
[35]), and the performance of several variations of DPCM has been tested by various 
researchers (Nijim et al., [54], Ross et al„ [63]). In this chapter we propose a 
compression strategy based on the finite differences for MR images. 



occurrences occurrences 
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For smooth images the k-th difference is sufficiently small and can be coded in lesser 
number of bits. The smallness of the k-th difference is used for compression. We 
quantitatively measure the smoothness of images by their membership in certain 
smoothness classes X ’ (defined more precisely in section 2 of this chapter). We choose 
k, based on the estimated smoothness of the image. Not much loss of information can be 
tolerated in medical image compression and these methods render lossless compression. 
In the present chapter, besides obtaining storage bounds for the proposed compression 
algorithm we present the error estimates for the representation of f. Rest of the chapter is 
organized as follows. 

We expect more compression through the finite difference methods, for images, which 
are sufficiently smooth. In section one, we discuss the method of storing k-th differences 
instead of pixel values for compression. Estimation of smoothness for the test images is 
carried out in section two. This is achieved by measuring the smoothness parameter a 
(defined more precisely in the following section). The compression algorithm and 
storage requirements for the proposed algorithm are presented in the same section. In 
section three a lemma due to S.N.Bemstein and certain other inequalities, which are 
useful in obtaining error estimates, are stated and proved. We present the direct error 
estimates and inverse theorems in section four. In section 5, the experimental 
implementation of the proposed compression algorithm is carried out on the test images. 
We show by experiments that for smooth images compression is more. 

3.1 Storage of Finite Differences: 

Throughout this chapter we will use the following notations. Let f be the digital image 
stored at a resolution m and depth b. Let N = 2 m and h = 2 m . The image f is stored as 
NxN matrix of pixel values, which are some averages of the analog intensity distribution 
function over small squares of size 2 m x2 m . Let f(m,n) be the mn-th pixel value. 
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Initial Storage: 

f(0,0) f(0,l) . . . f(0,N-l) 

f(l,0) f(l,l) . . . f(l,N -1) 

' ’ ’ ' ’ (3-1-1) 

f(N - 1,0) f(N-l,l) . . . f(N-l,N-l) 

For each fixed y, the image f(x,y) is a function of one variable. Throughout this chapter, 
we consider the image f as a function of y for each fixed x. Storage requirements for an 
image can be specified in terms of Bit rate, which is defined as follows 

Definition (3.1.1): 

Bit rate is the number of bits required for storing one pixel. 

Bit rate := No. of bits used for storage / No. of pixels used to represent the image. 

Bit rate is given as bits per pixel (bpp). The digital image of depth b has a bit rate of b- 
bpp. 

The k-th horizontal forward difference for the pixel data is given by, 

A° h f(x, y) : = f(x, y) 

and for j - 1, ,k A j h f(x, y) : = A j h ‘‘ f(x, y + h) - A j h _I f(x, y). (3.1.2) 

For a fixed k, the k-th horizontal forward differences, are calculated by using the above 
formula and the differences are stored as follows: 
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A” f(0,0) 

A' h f(0,0) . 

A^ f (0,0) 

Ah f (0,1) . . 

A k h f(0, N - k - 1) 

A'if(l,0) 

A' h f(l,0) . 

Ah f (1,0) 

A k h f(U) • . 

A k f(l, N - k - 1) 

A° h f(N - 1,0) 

A' h f(N - 1,0) . 

A k „f(N - 1,0) 

A k h f(N-l,l) . . 

A k ,f(N - 1, N - k - 1) 


(3.1.3) 


From these differences the actual pixel values for the image can be calculated using 

f ( r > s ) = X • AJ h f ( r >°) 0 < s < k - 1, 0<r <N-1 

j=0 \}J 

f(rJ) = A k h f(rJ-k) + ^](-l) u+1 [ k lf(r,j-o) k<j<N-l, 0<r<N-l 


U=1 




(3.1.4) 


The pixel values calculated by using the above expressions will be exact and there is no 
loss of information. We theoretically estimate the error between the original image f and 
the decompressed digital image. From the stored k-th differences, pixel values at non- 
nodal points can be estimated by using a k point equi-spaced interpolation scheme. 
These interpolation schemes help us to obtain images on a finer grid and are practically 
useful. At non-nodal points the following interpolation formula can be used 


f(x,y)=f(x,0) + ^A h f(x,0) + - 
h 


• + 


y(y- h )-[y-(k-i)h] 

k!h k 


A k f(x,0). 


(3.1.5) 


Here one argument in f(x,y) is fixed and yo,...,yk are nodal points. 


By image compression we are trying to reduce the redundancy in the original 
representation. For example, the Huffman coding uses the statistical redundancy, i.e. the 
number of occurrences of each symbol/state and can be calculated from the histogram of 
the data. The frequency information is already available with the data and this is used for 
variable length coding of the data. Our objective is to reduce the redundancy due to 
smoothness, which is a measure of correlation between adjacent pixels. With this 
objective we first estimate the smoothness, under suitable hypothesis on the images under 
consideration. Without any apriori knowledge about the smoothness of the images if we 
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store the k-th differences instead of the pixel values then it actually would take more 
space than the original representation. If the original image is coded in b bits per pixel 
then the k-th difference would require (b+k) bits per pixel. This is stated in the form of a 
lemma. 


Lemma 3.1.6. The k-th horizontal forward difference of a b-depth image takes (b+k)- 
bits at most. 

Proof by induction: 

1 difference: If we subtract two b-bit unsigned numbers (binary) then the result has at 
most b-bit magnitude and one sign bit which makes it (b+l)-bits. 

2 nd difference : If we take the seond difference using three b-bit unsigned numbers, then 
it introduces at most one carry bit extra and one sign bit. Therefore (b+2) bits in total. 

Let us assume that our claim is true for (r-l)th difference, the r-th difference 
A h f(x, y) : = A; -1 f(x, y + h) - A r h “' f(x, y) 
it is the subtraction of two r-1 th difference which already has a sign bit allocated and 
here it may at most introduce a carry bit extra, b+r-l+l=b+r bits are required at most. 

Although the k-th difference requires more space than the pixel values, the smoothness of 
the image allows us to code the k-th difference in less space than the original pixel 
values, for suitably chosen k. This is achieved by estimating the smoothness of images 
under suitable hypothesis. We fix the value of k based on the smoothness of the image 
and code the k-th difference instead of the pixel values. In the following section we 
estimate the smoothness of digital images by their membership in suitable smoothness 
class. The smoothness estimation gives the order of the k-th differences, this order 
information is used for compression. 

3.2 Estimation of Smoothness: 

We measure the smoothness of images by their membership in the following class of 


functions. 
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Definition 3.2.1 : For 0 < a < k, the smoothness space X a,k is the set of functions which 
have a smoothness, i.e. the k-th modulus of smoothness co k (S) is of the order of 8 a as 
8->0 

X“' k = {f :© k (f,8) = 0(6“)}. 

Let the image f e X“’ k , the value of a can be estimated by measuring co k (f,S), for various 
values of k. 

co k (f,8) = 0(8 a ), 


hence, co(f,S) < C 8 a , where C is a constant. Therefore, 

log(co k (f,8)) < log C + a log(8) . (3.2.1) 

The value of a and the associated constant can be estimated by measuring the slope and 
y-intercept of the graph log(S) Vs log(co k (f,5)). For estimating a we consider co k (f,5)i and 
w k (f,S )2 which are defined as follows 

• ® k (f»8), =max^|A k t f(i)| and 


oo k (f , §) 2 = max y|A k f(i)| 
k 2 i'i<yr' 1 


1 

2V 


For experimental estimation of a we take 8 = rh, r = 1,2,3 ,4, ... i.e. where h is the 
resolution of the given image, r is taken to be 8 steps. Graphs are drawn between 8 steps 
and co k (f , 8). The value of a and the associated constant C are estimated by measuring 
the slope and y-intercept of the graph between log(S) and log(co k (f,8)). The smoothness 
of the test images (test image set 1) are estimated. Throughout this chapter and the 
remainder of the thesis we consider images of depth 8 bits per pixel, i.e. images with 256 
gray levels. Plots of 8 vs coi(f,8)i and 8 vs ©i(f,S )2 are shown in Graphs 3.2. 1-3.2. 2. The 
plots between log(8) and log of ©i(f,8)i, ffl 2 (f,§)i, Wi(f,8) 2 , ©i(f,S) 2 are shown in Graphs 
3.23-3.2.6. For the test images the log plots are concave and the maximum slope is 
attained near the origin. The smoothness parameter with respect to co k for k=l,2 are the 



12 -omegal II -omegal 


maximum slopes attained in the corresponding log plots. The values of a and the 
associated constants are shown in Table 3.2.1-3.2.2. 


Graph of (5) Vs ® 1 (f,5)i 
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Graph 3.2.1 
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Graph of log(8) Vs Iog(oi(f,8)i) 
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Graph 3.2.3 


Graph of log(S) Vs log(cD 2 (fjS)i) 
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Graph 3.2.4 
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Graph of log(8) Vs log(o 1 (f,8) 2 ) 



Graph of log(S) Vs log(co 2 (f,5) 2 ) 



Graph 3.2.6 
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a-Estimates for test image set 1: 

Higher value of a indicates that the image is smoother and lower value of a indicates the 
image is less smooth. The estimated values of a for the test images are presented below. 
We consider functions which belong to X k ’ a , i.e. co k (f,h)=0(h a ). Along with the values 
of a, the associated constants C’s are also shown in the table. These constants refer to 
the constant involved in the O-notation. 


Images 

Smoothness measured in 11 

©i(f,5)i 

© 2 (f,S)i 

co 3 (f,6)i 

co 4 (f,6)t 

a 

C 

a 

c 

a 

C 

a 

C 

60.img 

0.544 

1322.5 

0.586 

1709.0 

0.529 

2902.0 

0.490 

5292.1 

64.img 

0.559 

1309.8 

0.664 

1703.5 

0.664 

2847.6 

0.654 

5137.7 

68.img 

0.555 

1442.0 

0.658 

1901.4 

0.669 

3182.2 

0.660 

5751.3 

72.img 

0.556 

1653.7 

0.646 

2178.2 

0.648 

3658.1 

0.651 

6610.5 


'able 3.2.1 


Images 

Smoothness measured in 12 

coi(f,5) 2 

©2(f,8)2 

G> 3 (f,§)2 

co 4 (f,5) 2 

a 

C 

a 

C 

a 

C 

a 

C 

60.img 

0.743 

9.541 

.0.942 

10.45 

0.906 

16.66 

0.855 

29.63 

64.img 




11.77 


18.23 



68.img 

0.717 

11.38 

0.972 

12.96 

1.038 

20.27 

1.051 

20.27 

72.img 

0.666 

11.99 

0.832 

14.59 

0.871 

23.67 

0.889 

41.84 


Table 3.2.2 


The value of oc gives an indication about the number of preceding zeros in the k-th 
differences. These preceding zeros need not be stored and consequently the memory 
space requirements for storing the k-th difference is less when compared to storing the 
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actual pixel values. Here the compression is done precisely by not storing the preceding 
zeros. At the lime of decompression, the numbers of preceding zeros are known from the 
value of a. Now we present the compression algorithm. 

Algorithm 3.2.1: 

Step 1 . Estimate a and the associated constant for the given image using the relation 
(3.2.1). 

Step 2: Choose k such that 0 < b+k+l-ma < b. 

Step 3. Calculate the horizontal forward differences of the image recursively by (3.1.2) 
for the chosen k. 

Step 4: Store differences available in the form (3.1.3) using Huffman coding (cf: 1.2.1). 
Storage requirements: 

Let f € X k,a . The first k differences for each row are stored separately. If f is stored in b 
bits per pixel then each k-th difference A k h f requires b+k bits for storage. Out of b+k bits 
one bit is reserved for storing the sign. We have, 

A l l f(x)|£ffl l (f,h) = 0(h-). 

For h = 2' m the resolution of the image, we have |A k h f(x)| = 0(2“ ma ) , therefore the binary 
representation of the k-th differences will have ma preceding zeros which need not be 
stored, therefore the net bit rate required for storing the k-th differences is b+k-ma. 

The initial (k-1 ) differences occurring in each row of (3.1.3) are stored in the header and 
the rest of the k-th differences are coded in b+k-ma bits. Compression is achieved 
because (b+k+l-ma) < b, this follows from the way k is chosen (as per the algorithm 
3.2.1) based on the estimated a. Consequently the total storage space required for NxN 
size image will be (N-k) x N x (b+k+l-ma). Compression through algorithm 3.2.1 is a 
non-lossy. At the nodal points we can exactly get back the original pixel values using 
equations (3.2.1). The exact implementation of algorithm 3.2.1 is discussed in section 


3.5. 
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3.3 Bernstein’s Lemma and certain estimates for smooth images: 

Here we estimate the error incurred in k-th difference decompression by measuring the 
difference between the original, analog intensity distribution function, and the 
decompressed digital image. We started with a digital image, which are samples of the 
original function. After decompression the values at any non-nodal point is estimated 
using a local k-th degree Lagrange interpolation polynomial (see Hildebrand [28]), L k . In 
this section we shall prove some lemmas which are used in obtaining error estimates. 

Lemma 3.3.1 (S.N.Bernstein) (seeTiman [71]): 

If the continuous function f(x) on the interval [a,b] has (k+1) distinct zeros, then its k-th 
difference A k f( x) has at least one zero on this segment for sufficiently small 8. 

Since the original reference of Bernstein was not available to us an inductive proof of the 
lemma is included below: 

Between any two adjacent zeros of f viz. z k zj we will show that there exists a ^ e (zj, 
Zj+i) such that Ag. f (£,) = 0 for sufficiently small 8j, 

%*■,«)> 0 (<*< 0 ) 
because Zj and Zj+i are adjacent 

f I, , attains local maxima (or local minima) say M. Consider lines y = r (0 < r < M). 
For each r, y = r intersects with f|( ZjiZj>I ) at least at two points except possibly for the case 
r = M. 

Let pj and p j+! be two adjacent points of intersection of y = r and f| (ZjjZj+l) • Without loss 
of generality let pj+i > pj so that 0 < 5j = pj+i - Pj • For this 5j , A s . f(pj) - 0 and V 8 < 8j 
we have Ag. f(p) = 0 for some p e [ z j> z j+i]- 

=> Between any two zeros Zj , Zj+i of f there exists a zero of A s . f(x) . 
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Let 5i - min{5j} 0 < j < k. For this Si 3 at least one point £j between every pair of 
adjacent zeros of f 3 Ag. f(£, j ) = 0 zj < Zj+1 . 

=> A s . f has at least k-zeros. A' 5 . f(x) is a continuous function because 
Ag f(x) = f(x + 5) - f(x) and f is a continuous function. 

— ■ > ^ * s a continuous function and has k-zeros. Therefore by the same above argument 

3 a 5, >0 such that between any two zeros of A' 5 .f there exists a zero of A' g , (AV f ) . 
For 82 = min {5i, 5 2 }, A' 8j (A' 5j f(x)) = A 2 5j f(x) has a zero between any two zeros of 
Ag. f continuing this process we get, there exists a £ o A k 5 f(£) = 0 where 8 = min {8j} 0 < 
j < k i.e. 8 < min {| Zj - Zj+i |} = minimum distance between any two adjacent zeros 

'Fhe following lemma gives a measure of the remainder term in interpolation by n-th 
degree polynomial in terms of (n+l)th differences. 


Lemma 3.3.2: 

For any real continuous function f(x), the following Newton’s formula is valid 


f(x) = f( 0 ) + — A h f( 0 )+- 


•+ X(X h> n ^„ — " 1 ~-A n h f( 0 )+R,.i (x) . 


where R,.,(x) = with 0 < 8 <h. 


Proof : 

Let, P(x) = f(0)+^A h f(0)+ 


x(x - h)- • • [x - (n - l)h] 
+ n! h” 


A h f (0) 


Consider g(t) = f(t)-P(t)-[f(x)-P(x)]^_ x °| (x-x n ) (3.3.1) 

Equation (3.3.1) is continuous and has (n+2) zeros viz., t = x 0 ...xn, x, where the points 
Xo...x„ are equi spaced (spacing h) and x is at most h distance away from one the the x/s. 
Therefore the maximum distance between any two adjacent zeros of g(t) is h. By Lemma 
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3.3.1, for the function g, 3 a £ 3 A" 5 +1 'g(£) = 0 and 0 < 6 < h. Applying A 1 
the point £, we get 

o= AV'gtO - AVf(4)-[f(x)-p(x)] 7 ( n+1)!5 ‘" hence 

l X X o)**’( X X n ) 

f(x> - 


" +1 to (3.3.1) at 


(n + 1)!S 


x(x- h)---[x-(n-l)hl , 

^“ f © 


[0 < 5 < h]. 


Now we present the interpolation process, which is adopted for evaluating the pixel 
values at non-nodal points. 

Process: 

Here we process each row of the image separately. The function under consideration is 
that of one variable, keeping the other one fixed and here interpolation refers to 
construction of interpolation polynomial using equi distant data. Let f e C 2 % and be 
approximated by k-th degree algebraic interpolation polynomial L f k . At finitely many 
equi spaced points the value of f is known. We partition I=[0,2 tc] into Ii,...,I n such that 
each Ij contains (k+1) nodes. In each Ij, f is approximated by L kj , the k-th degree 
interpolation polynomial constructed using the (k+1) equi spaced nodes. 

L f k,h (x) = IX 00 L k,hj( x ) ( 3 - 3 - 2 ) 

i=l . 

L W x ) = ZW x ) f ( x J+> 

i=0 

(x — Xj +0 )‘ **(x — x j+i _ 1 )(x - Xj +i+1 >-(X" Xj+k) ^ 2 2} 

lj+i(X) = (X j+i -x j+0 )-(x j+i -x j+M )(x j+i -x j+i+1 )-0t j+i -x j+k ) ' 

L f k 2 _ m is the interpolating polynomial on the whole of f. X,. is the characteristic function 
on Ij. This amounts to piece wise approximation of f by interpolation. 
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Theorem 3.3.1: 



Jf L' k h is the decompression of f e [x 0 _h,x k +h] 
(h = 2' m ) is used for compression then there exists a 


, whose digital image at resolution 
S, 0 < 5 < h such that 


hi 



Proof: Follows from Lemma (3.3.2). 



Remark: 

The right hand side in the expression (3.3.4) has This could be large f 0r 

general functions, however for the images under consideration (test image set I) we Se 
from the following graphs 3.3.1-3.3.2, that the slope is maximum near the origin and t he 
shape of the graph is concave. Therefore for the test images under consideration a bo Ul J 
based on the slopes could be obtained from graph of e> k+1 (f,8) vs 6 k+1 . For k=l,2 th e 
value of right hand side comes out to be approximately 0.03. Theorem 3.3.1 says that f 0r 
k = 1,2 the error incurred in the process is around 3 % for the test images und er 
consideration. In the next section, we obtain the precise order of error in o Ur 
compression scheme. 



12-omega 2 


Graph of 5 2 vs co 2 (f;8) for test images 



delta square 


Graph 3.3.1 


Graph of 8 3 vs © 3 (f;5) for test images 



Graph 3.3.2 
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Lemma 3.3.3: 

If f is a bounded function in [x 0 -h,x, +h] then the k-th degree Lagrange interpolations, 
L k , are uniformly hounded operators, independently of h. 

Prof: 

L U X ) = Z 1 - ( x)f(x 1 ), 

i - 0 


where l ; (x) = 


(x-x 0 )---(x-x i _ 1 )(x-x i+1 )-"(x-x k ) 

(x i -x 0 )...(x i -x w )(x i -x i+ 1 )-(x i -x k )' 


C h k 

Hence, |lj(x)| < — < C 3 
v -> n 



where C k ,is a constant independent of h and depends only on k. The result follows from 
this. 


3.4 Direct and Inverse Theorems: 


In this section we derive direct and inverse theorems for the process presented in the 
previous section. The direct theorem is derived using Peetre’s K-functional, its 
equivalence with the modulus of smoothness, and certain properties of modulus of 
smoothness whose short proof for the Lp-version, adapted from the sup norm case in 
Lorentz [46], are included for the sake of completeness. 


Definition (3.4.1) (Peetre’s K-functional): For 1 < p < °o, the Peetre’s K-functional is 


the quantity 





This is equivalent to the usual (k+l)-th modulus of smoothness C0k+i(f,t)p. This is stated 
and proved as a theorem. 
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are independent of t 


Theorem 3.4.1 : For 1 < p < oo the Peetre’c v fi lr w , • 

’ Cetre S K - fu nctional 1S equivalent to the modulus of 

smoothness, K( f , f).. ~ f t \ : _ 3 „ 

v , vp <»ku( r , tj p i.e., 3 constants C u C 2 > 0 and 

such that 

Ci co k+1 ( f,t) p < K(f,t) p ^C 2 co k+1 (f,t) p . 

To prove this we need the following lemmas, 

Lemma 3.4.1 : For 1 < p < oo, the inequality © k (f,t) p < 2 k ||f|| p j s valid. 

Proof: For 1 < p < oo, 


2 3X 


o k (f,t) =sup flA k h f(x) P dx 

N<t 


= sup| 

1h|<t 


2tc ^ 

jE(-i) 1 -' 

0 r=0 




v r y 


f(x +rh) 


dx 


; E 

r=0 V 


f 2« 


f(x + rh) | p dx 

Vo j 


V 


-2 k Ilf 1L 


For p = oo, we refer Lorentz [46], 


Lemma 3.4.2: For 1 < p < oo, f e C k and f (k) e Lp[0,27i] the inequality 
<t k |jf (k) || p is valid. 

Proof: For 1 < p < oo, 

l t 

J-- Jf <k) (x + y,+-”+y k )dy,---dy k 


0 0 


V o 


A k f(x)| P dx = J {••• /f <k) (x + yi+-"+y k )dy,-"dy 1( 


0 0 


dx 


t 


\ 

< !••• (x + y 1 +---+y k | P dx 

0 0^0 ' 



dy, •••dy k (By Minkowski’s ineq.) 


<t‘ 


We now present the proof of Theorem 3.4.1. 
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Proof of Theorem 3A1 : 

Let f a „ be the (k+l)-th Sleklov mean (its definition and properties are presented in 
chapter 2). 


K(f,t) p < 



t.k a i 


+ t k4 ' 


If <k + I) 
t,k + l 


£ + t k %t-< k+ x +1 (f,t) p 

* C 2 co k4 i (f,t) p 


[Where C 2 = 2 max {1,C,}] 


Conversely, we choose g <= C k+1 
co k+1 (f,t) p < co k + 1 (f - g, t) p + co k+1 (g, t) p 

- ^ Ik - dip + t k+1 ||g (k+1) | p [By lemma 3.4.1 and 3.4.2] 

< 2 l,l {l|f-g|| p +t l *'||g <k *' , || p } <2“'K(f,t),. 


Remark: 

We prove the direct and inverse theorem for a general function for the case p = co. For 
the case 1 < p < oo, the Lagrange’s interpolation operator is not a bounded operator over 
C 2w . However, an MRI image corresponds to a trigonometric polynomial in two 
variables. Our present compression implementation corresponds to the rows of the image 
matrix giving rise to trigonometric polynomials in one variable, namely x, of a fixed 
order. Since the space, 3„ of trigonometric polynomials of order n is finite dimensional 
for the MRI purposes, the point functionals T n (Xk) are Lp-bounded: 
|T„(x k )|<|T„*D n (x k )| 


[ D n stands is the Dirchlet kernel and * stands for convolution.] 



with 1 < p < co and — + — = 1 

P P 


[By Holders inequaltiy] 
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Let f be a MR] image, then for I < p < oo the k th r 

T , , , . ’ k th degree Lagrange interpolations, L k , are 

Lp-boundcd opciators, independently of h 

Proof: 

By .be remark above .bis ,e mma> we see ,ba< the pen, toctlonals for ^ ^ purposes 
are L.-boam.ak ! be proof of ,be above lemma f„ 1Iows on the same ]ines of ^ 


Lemma 3.4.4: 

g C [\ () h,x k +h]. If L g kh be the k-th degree Lagrange interpolation 
polynomial of g, then for 1 < p < oo, the following inequality holds: 


y-f-u, ^ C k h k+I g (k+1 > 


where C k is a constant depending only on k. 


Proof: 


The result for the case p = oo, i s well known (see Hildebrand [28]). We shall prove the 
result for p - 1 . From these two cases, the result for all p, (1 < p < oo) follows by the 
apphcation of the Riesz-Thorin interpolation theorem (see Bergh and Lofstorm [5]). 

Let P k be the k-th degree Taylor’s polynomial of g, about a point a e [x 0 - h, x k + h], say. 
Then we have, 

||s- l '4 HM -Ll.il, Sfc-P.l dlLSl, 

Using the integral form of remainder in the Taylor’s theorem, we have 

***** 1 x k +h X 

«-M. = Jlg(x)-P k (x)|dxsi { I(x-s)|*|g ( ‘* l, (s)dsdx 


x 0 -h x 0 -h 
x t +h 


Xt +h 


< 


J |g <k+1) (s)| |(x-s)| k dxds 


x 0 -h 


Xo-h 


k! 

<i((k + 2)h) , * l |g l ‘* l »( | 


<C t h“|g 


k+I L(k+1) 


The constant C k depends only on k. 
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Theorem 3.4.2 (Direct): 

If L' Ul is the decompression off e C 2lt , at any arbitrary point of f whose digital MR image 
at resolution m, re. h 2 is used for compression then the following estimate holds 


f-U 


k,h 


^Cco k+I (f,h) D , 


where the constant ( is independent of f and is dependent only on k. 

Proof: 

We consider a function g e C^' with g 0 ^ continuous and j|g (k+1 >|| p is finite , we have 


f-U 


k.h 


f - g + g L k h + L k h L kh 


£ |f - g + g - L g k 


"k,h 


+ 


r g - f 

H.h 


4 +»||)l f -4 

- C t K(f, h) p 
S0> lfl (f,h) p 


+ C,h 


k+I 


jCk+1) 


The constant C\ depends only on k. 


[By lemma 3.4.4] 


[By theorem 3.4.1] 


Theorem 3.4.3 (Inverse): 

For 1 < p < oo, the following inequality holds 

®k + i(Uh) p <C 

where the constant 0, is independent of h and f. 


f-U 


k,h 


Proof: 

For a given h > 0, we construct Lagrange’s interpolation polynomial L f kh from the 
available samples of f. L k h is an algebraic polynomial of degree k, therefore 
» k .,(f,h), = © t ., (f - L k h , h) p < 2 l+1 |jf - L k h || p . 

Combining the results of Theorems 3.4.2-3.4.3, it follows that the error in our 
interpolation scheme is precisely of the order of © k+ , (f,h) p . 
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3.5 Experimental Results: 

Algorithm 3.2.1 is implemented on the test images and here we present the results after 
loss less compression. Algorithm 3.2.1 is implemented in integer arithmetic as follows: 
We consider test images at depth 8 bits per pixels. Image data is available in a 256 x 256 
grid. Pixel values are available at spacing h =1/256. Let the range of pixel values in the 
test images be normalized to take values from 0 to 2 8 -l. Let {Py} be the pixel data. The 
horizontal k-th differences for suitable k are calculated for the pixel data using the 

equation (3.1.2). After computing the matrix of differences (cf: 3.1.3) corresponding to 
the pixel data they are stored as follows: 

(i) The entries of difference matrix are signed integers. Let the entries of the 
difference matrix be denoted by {d,j}. Let r be the smallest negative number in 

the difference matrix. 

(ii) dy = d,,+ j r | , are computed for the difference matrix. This is done to avoid the 

storage of a sign bit, required for each entry in the difference matrix. 

(iii) The storage matrix along with the header containing the value of r is entropy 
encoded using Huffman variable length coding (cf: 1.2.2) 

The resulting images after decompression are shown in Image 3.5.1. We note that the 
entropy of the storage matrix djj is less than the entropy of the actual image. Along with 

each of the decompressed image, the corresponding percentage compression, bit rate 
required for storage and the entropy of differences is given. We note that the performance 
of lossless compression algorithm 3.5.1 is comparable with JPEG coder at minimum loss 
levels (Image 3.5.2). The results of algorithm 3.5.1 applied on sharpened images are 
shown in Image 3.5.3 and Table 3.5.1. The smoothness estimates for the sharpened 
images are also shown in Table 3.5.1. 



N on-lossy Compression through Algorithm 3.2.1 
Images after decompression 



60,img decompressed 
P.C.: 40.30 % B.R : 4.77 bpp 
Entropy of differences : 4.682 



64.img decompressed 
P.C.: 41.08 B.R.: 4.71 bpp 


Entropy of differences : 4.597 



68.img decompressed 
P.C.: 38.83 % B.R.: 4.89 bpp 

Entropy of differences : 4.757 



72.img uncompressed 
P.C.: 36.03 B.R.: 5.11 bpp 

Entropy of differences : 4.991 


Image 3.5.1 






Image compression through JPEG 



60.img.jpg 
Bit rate = 4.67 bpp 


64.img.jpg 
Bit rate = 4.66 bpp 



68.img.jpg 
Bit rate = 4.91 bpp 


Image 3.5.2 


72.img.jpg 
Bit rate = 5.13 bpp 






Compression of Sharpened images 



img.36 


img.3612 



img.36 123 



img.36 1234 


Image 3.5.3 


The percentage compression and bit rates are presented in the following table. 
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in (3.1.2). 
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From the table ' .1 1 -o,2.2 we note that the smoothness of the test images decreases from 
img.60 to ime 72. And we note that the compression is more for images, which are 
smoother. The main advantage of algorithm 3.5.1 is that it is non-lossy which is highly 
useful in medical imaging. In chapter 4 of the present thesis, we develop an algorithm 
based on reversible smoothing, which work well for the images under consideration 



CHAPTER IV: COMPRESSION THROUGH 
RF.VERSIBLE SMOOTHING 


4.0 Introduction: 

We hav e seen in the previous chapter that for smooth images the k-th differences are 
sufficiently small and can be used for compression. In the present chapter we propose 
and analyze a compression algorithm by reversible smoothing. We obtain a smooth 
version of the original image by convolution and store the smooth image by the k-th 
difference encoding as discussed in the previous chapter. At the time of reconstruction 
we first obtain the smooth image from its differences and an approximation to the 
original image is reconstructed from the smooth image by deconvolution. Let a given 
MR1 image be, 

f(x,y) = £c. k e i(jx+ky) . (4.0.1) 

We apply an invertible transformation on the Fourier coefficients to get a new function 
? . This transformation could be based on a known convolving function or otherwise. 

For the chosen convolving function, 

7 (x, y) = £ c jjk e ,(J x+k y) . ( 4 - 0 - 2 ) 

j.k 

f is the original image and 7 is the smooth version of f obtained by convolving f with 
suitable approximation kernel. In practice Cj k ’s are obtained by multiplying Cjx s with 

p(k)’s, the filter values corresponding to the convolving function. Instead of storing f by 
algorithm (3.2. 1 ) we store 7 by the same algorithm and at the time of decompression we 
get 7 first and from that f can be obtained by deconvolution. The process of 
compression through reversible smoothing is illustrated using the following example. 





Original - Image 4.0. 1 Smooth - Image 4.0.2 

P.C: 34.84% B.R.: 5.21 b P.C: 48.37% B.R.:4.13 


Decompressed: Image 4.0.3 
Relative 11 error : 0.0618 
Realtive 12 error : 0.0472 
Average 11 error : 0.0074 
Average 12 error : 0.0092 


Error Image 4.0.4 

Absolute difference between Image 
4.0.1 and Image 4.0.3. Intensity 
increased 25 times 


P.C - Percentage Compression, B.R-Bit Rate. 
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We consider i'.»ur different approximation processes Fejer, Poison, Jackson and Steklov 
means for smoothing the original images. The rest of the chapter is org ani zed as follows, 
the four well known approximation processess listed above are reviewed and the 
corresponding filter coefficients are presented in the first section. In section two, we 
discuss the compression strategy and study the effect of smoothness in compression. 
Experimental implementations of the compression strategy are also carried out in section 
two. In section three a relation connecting sampling resolution and the correspon din g 
depth is obtained. Section four is devoted to wavelet based compression methods. 


4.1 Approximation Processes: 

Let the original image f e C:* , and let it be represented by n-th partial sum of Fourier 

series, 

S n (f;x) = ~ + X a k coskx + b k sinkx (4.1.1) 

2 kssI 

where a k - ■■ jY(x) cos kxdx and- b k =— |f(x) sin kxdx. 

^ ft ~7t 

It is well known that S n (f) does not converge to f even for continuous functions. Several 
approximation processes are available in the literature (Zygmund, [79], Natanson [53]). 

Most of them are available in the form of ~ + Xp(k)( a k C0S k x + K sinkx)* where 

■4 k=l 

p(k)’ s are the corresponding filter coefficients. A brief description of the filter 
coefficients corresponding to the approximation processes by using Fejer, Poisson, 
Steklov and Jackson kernels are as follows. 


Fejer: 

These are the arithmetic averages of the n-th partial sum 

■ 2 nt 
sm — 

obtained by convolving f with the Fejer’ s kernel ~ ■ 

sin 2 - 


of the Fourier series and are 


F n (f) = 


S!+Sf +•••<_ 


n 
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, . S1 „=5^ 

F„(f;x) = - — J — 2 — 

2 


, . . . ^ jj \ ' n k | / , 

! 'n (,;x > = 2 + M n ~JlA coskx + b k sinkx) 


(4.1.2) 


The factor 


n k \ 


acts as a weight on high frequency terms and kills the high frequency 


components thereby introducing some amount of smoothness in the image. 


Poisson: 

Poisson sums are obtained by convolving f with the poisson kernel. For 0 < r < 1, the 

Poisson kernel P(r, t) is given by 


P(r,t) ~ 1 + 2]>V cosjt = 


l-2rcost + r 


B n ^ 

P„ ( 0 * + X r k ( a k cos kx + b k sin kx) 


J 

** k»l 


(4.1.3) 


Here, the filter coefficients are generated by using r k (0 < r < 1). Smaller the value of r 
higher the smoothness, practically r cannot be too small because of the roundoffs that 
occur at the time of deconvolution, consequently we may loose some information, if the 
image is smoothened too much for the sake of compression. We show by experiments 
that there exists a practical choice of r, which yields good compression and acceptable 


roundoffs. 

Jackson: 


Jackson's sums are 


» obtained by convolving f with the Jackson kernel given by, 


1 sin o 


2?xn . t 
sin 


( n(t - x)^ 
it sin 

, , r . . _ -J f ^ f( t )dt 

J,,(f ’ } ~ 2rcn(2n 2 + 1) .[ sin l^ 

A 2 * 


(4.1.4) 
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Rathorc |M] has t’lu'n an explicit construction of a sequence of linear trigonometric 
polynomial operators of the host approximation order as given by Jackson’s theorems 
(see Natansonj \t] 1 \ particular ease (Lemma 4.1.1) ot his general theorem can be used 

to obtain the filter coefficients in the case of Jackson operators. 


Lemma 4.1.1 : 

The following identity holds true 


nt 


sm 


sin 


L+ £l k coskx, 


where L 4n' 


1 Ldf, k 

4 > 1 - “ 

2 f-fk n 


4n + 2n 


k •. 


and 


„ , , k(K l)(k 1 1) 2(n - k)(n - k - l)(2n + k - 1)) 

2n(k i 1) t , + — — — 

, I 3 3 

k = j (2n - k)(2n k -l)(2n ~k + l) 

3 


1 <k <n 
n < k < 2n - 2 


The filter coefficient p(k)’s for the Jackson sums are obtained as follows, 




t 

2n' t 1) J 


2nn(2n' t 1) 


. n(t-x)V 

sm 

2 

t - x 

v sin “r ^ 


f(t)dt 


2itn(2n‘ + 1] 


n 2n«2 

L | !ft)dH ft 2n(2n 2 + 1) 


it % 

Jf(t)coskt coskxdt+ Jf(t) sinkt sinkx dx 




31 1 


n(2n 2 +1) 2 2n(2n 2 +1) 

where L and L arc as in the Lemma 4.1.1. 


ja k coskx + b t sinkxj, 


(4.1.5) 
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Steklov: 

Ther-th Steklov t I rniun |71] ) means of f at a spacing h, f h>r are given by 
f ht (x) * f 1'*, , i x ) dx with f h ,o(x) = f(x) 


(4.1.6) 


Lemma 4.1.2; 

If (4.1.1 ) is the partial sum of the Fourier series of f, then for any integer r > 1, 

-A, ( si n khV 


f, . . (x) 


o I 11 Ml / \ 

4 L v — kh — J cos kx + bk sin kx ) (4L7) 


Proof: 


Let r ■ 1 . f, ( \ t 


% * h 


f(x) dx 


2h,C'"'”" 2h ^ 2 k „ 


i k 4 h. 

1 r( a 
2h 


" + X ( a k cos kx + b k sin kx) dx 


x-ii 


n j j oh J X*b 

i V a k “ |coskxdx+b k — Isinkxdx 

fcl { 2h , J h 21 l x-h 


x 

c 

J 


a„ A ( s in kh V . . . , i 

" } \ [a k coskx + b k smkxj. 

k»i V kh J 

For r = 1 , (4.1.7) is true let us assume that it is true r = m-1 and let us prove (4.1.7) for 


r = m 


Mi. m 


(X) 


2h 


r " 1 x+ fY sin kh 

J Vm |1 X ) ~ Oh •* V 


2h s V kh 


m-1 


^ n 

” + X ( a k cos kx + b k si n 
V 2 k=i 


2 ft* kh ) 


n / . • , \ n* If , X4h 1 x+ ^ 

a » „ W Sll| kh \ I _ J coskxdx + b k — Jsinkxdx 


\ 


2h 


x-h 


x-h 


*“ i (a k coskx + b k sinkx). 


2 f-fi kh 


Hence 4.1,7 is true by induction on r. 



4.2 Compression Algorithm: 


We now present the compression algorithm based on reversible smoothing, let the given 

MR1 image be of the hum (4.0,1 ), 

4.2.1 Compression: 

Step l:Fix the approximation process and calculate the filter coefficients p(k)’s. 

Step 2: Compute c yV s in (4,0.2) using the calculated p(k)’s and c j)k ’s and obtain smooth 

image f . 

Step 3: Apply Algorithm (3 2, 1 ) on f , and entropy encode the corresponding k-th 

differences using Huffman’s variable length coding (1.2.1). 

4.2.2 Decompression: 

Step 1 : Decode the variable length coded file by the procedure described in (1.2.1) and 
obtain the k-th differences of the smooth image. 

Step 2:0btain the smooth image from its stored k-th differences, using the relation 

(3.1.4). 

Step 3: Reconstruct the image from the smooth image by deconvolution by the same 

convolving function used for smoothing. 


The filter coefficients (windows) corresponding to the approximation processes listed 
above are as shown in the Graph 4.2.1 . For the case of Poisson the value of r is taken to 



WINDOWS 



Experimental setup: 

Here we smooth the 256x256 test images using the processes considered above. The test 
images have 256 gray scales, which take 8 bits per pixel for the original storage. The two 
measures we use for estimating the amount of compression are percentage compression 
and bit rate. Here the initial (original) storage for a 256 x 256 image at a bit rate 8 bpp is 
65536 bytes. Although the transformation applied on the Fourier coefficient is invertible, 
by dividing the filtered coefficients with the respective filter values, it leads to some 
round off errors. So there is a practical range beyond which images cannot be 
smoothened for the sake of reversibility. 

In case of Poisson approximation the value of r determines the imposed smoothness. 
Lower the value of r higher the smoothness obtained, but r cannot be too low. To 
estimate the useful range of r for the class of images under consideration, we measure the 
error after decompression. Four different quantities are used (cf: 2.6), relative 11 (rll), 
relative 12 (rl2), average II (all) and average 12 (al2). Besides these error measures to 
analyze the performance of the compression algorithm we introduce a measure of gain, it 
is defined as follows. 
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ill cam 


Compression Ratio 
relative 11 error after decompression 


Compression Ratio 
relative 12 error after decompression 


(4.2.1) 

(4.2.2) 


As r decreases (torn 1 .*). the compression ratio increases as the smoothness of the image 
increases I he lelaine enoi decreases up to some point and then increases, the point 

where relative gam attains maximum is the region of interest. 

Result of the compression algorithm 4.2.1 using convolution through Poisson kernel is 
shown in Graph 4.2 .2 -4.2 \ Variations of compression with respect to smoothness in 
terms of relative gam arc shown in Graph 4.2.6-4.2.9. It is seen that for the convolution 
of pixel data using Poisson kernel, the value of r = 0.985 works well for the test images 
under consideration Images tiller decompression for the case of Poisson, Fejer and 
Stcklov are shown m Image 4.2.1 -4.2.3 and 4.2.5. Sharpened images (Test image set 2) 
arc compressed using the algorithm 4.2.1 and the results are presented in image 4.2.4. 
Results indicate that for visual usage of images compression strategy 4.2.1 based on 
Poisson process is better suited. The compression ratio attained for the test images are 
higher when compared to other methods studied in the previous chapters of this thesis 
















V ariation of gain with respect to smoothness (Poisson process) 




Graph 4.2.7 




Graph 4.2.9 


In case of Poisson approximation, the graph of Smoothness(r) Vs gain indicates that there 
exists a point at which the gain is maximal The round off error at that value of r is 
admissible (there is no visual loss of information in the images obtained after 
decompression). As r decreases, the gain increases initially but after some stage the 
round off errors in deconvolution dominates and hence the gain starts decreasing 





Decompressed images corresponding to maximal gains for the test images are shown in 
Image 4.2.1 in the case of Fejer approximation, the k-th Fourier coefficient f(k) is 

N- j {£ j 

multiplied by a quantity — (filter coelficients) where N is the highest frequency, 

(in our case it is 128) and k runs from —127 to 128. During the process of 
decompression, the Fourier coefficients calculated from the smooth image are divided by 
the filter coefficients. When the value of k is 128, we see that the multiplier becomes 
zero and causes problems at the time of decompression. To avoid this, we take N=129 
and compute the filter coefficients. We observe that for N = 129, the amount of 
smoothness imposed by filtering is more. Although, the smooth image obtained by 
using, N=129 gives better compression, it lacks reversibility. The decompressed image 
4.2.2 (corresponding to N=129) in this process shows up deterioration. As the value of N 
increases, the visual deterioration starts decreasing. For all practical purposes, we use the 
value of N 130, which strikes a good balance between compression and reversibility. 
The graph of N vs decompression error is shown in Graph 4.2.10-4.2.11 and the images 
after decompression for various values of N are presented in Image 4.2.2-4.23. In case of 
the Steklov approximation, not much smoothing occurs. The Steklov window as shown 
in Graph-4.2.3 is almost horizontal and all the filter values are near 1. The decompressed 
images are comparable to that of images obtained in Poisson’s process. 

In comparison with Fejers and Steklov, the Poisson approximation process works better 
for the test images considered here. The additional advantage of using the Poisson kernel 
is that the value of r can be suitably adjusted, so for different values of r we get different 
compression. We present the images after decompression through these processes along 
with their compression percentage and the respective errors. 
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Poisson process - Images after decompression 




Percentage compression : 69.94 
Bit rate : 2.40 bpp Size : 19697 bytes 


Percentage compression : 66.66 
Bit rate : 2.67 bpp Size : 21852 bytes 



Percentage compression : 66.46 
Bit rate : 2.68 bpp Size : 21982 bytes 


Percentage compression : 59.33 
Bit rate : 3.25 bpp Size : 26656 bytes 



Image 4.2.1 






average II -error 


ny 


Variation of Smoothness with N in Fejer’s smoothing 



Graph 4.2.11 
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Fejer smoothing - Images after decompression for various values of N 




N=130 N=131 


Image 4.2.2 






Image 4.2.3 








COMPRESSION OF SHARPENED IMAGES (Using test image set • 21 
POISSON PROCESS-ALGORITHM 4.2 g l § ' 2) 


Images after decompression 



r = 0.9846 P.C : 63.94 % B.R.: 2.88 r = 0.9837 P.C : 61.61 % B.R.: 3.07 
Relative error rl 1 : 0.0567 Relative error rll : 0.0628 



r = 0.9828 P.C : 60.62 % B.R.: 3.15 r = 0.9819 P.C : 60.53 % B.R.: 3.16 
Relative error rl 1 : 0.0749 Relative error rll : 0.08 1 1 


n the above images, r is the corresponding Poisson kernel parameter. 

Image 4.2.4 






Steklov Smoothing 



Steklov smooth Image corresponds to 
img 60 Compression : 67.04% and Bit 
rate = 2.64 bpp 



Decompressed from Steklov smooth image 
Relative 11 error =0.076: Relative 12 
error = 0.063: Average 11 error = 0.00095: 
Average 12 error = 0.013 


Image 4.2.5 

In the next section we study the effect of smoothness in sampling. The resolution and 
depth of a digital image are related to the smoothness of the image. 

4.3 Non Redundant Sampling: 

The depth and the resolution of a digital image determine the initial storage requirements. 
If the image is sufficiently smooth, then for a required depth, the resolution can be coarse. 
The extreme example is a constant image which requires just one pixel information to 
reconstruct the entire image. Detailed images require fine resolution. The quantity 
modulus of smoothness plays an important role in fixing the resolution for a given depth. 
For a given depth, the resolution and smoothness measured in terms of modulus of 
smoothness can be related. 
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Definition 4.3.1 : I he Space X K *■“ is the set of functions whose co k+1 (f,h) is of the order 

ofh", where 0 < a < k+1, 

X Mi ' a = { f : co k+ i(f,h) = 0(h a ) }. 

Let I 1 1, be the k-th degree inteipolution polynomial (cf: 3.3.8) of f, the analog intensity 
distribution function. L' u , is constructed using the k-th differences stored as per the 
compression algorithm (3.2.9). L kh is said to be a b-perfect decompression of f, if the 
first b-bits of L f k h matche with that of f. Resolution m (h = 2' m ), and depth b bits per 

pixel are two main parameters of a digital image. In this section, we address the question 
of what is the minimum resolution required to reconstruct the image at any arbitrary point 
in the domain under consideration, such that the reconstructed image is b-perfect. 
Analogously, the question of what is the maximum possible depth that can be attained at 
any arbitrary point in the domain of the image, by using a digital image at a resolution m. 
The pixel value at any arbitrary point is evaluated using L k h (with .h = 2' m ). The 

following theorem gives the relation between depth and resolution of a digital image in 

terms of the modeled smoothness. 

Definition 4.3.2 (Maximum attainable depth b m using L f k 2 _ m ): 

Let f be the original image and L f k , „ be a reconstruction of f from the pixel data at a 

resolution m (spacing - 2“ m ), by using a local Lagrange’s interpolation polynomial of 
degree k. Then, the maximum attainable depth is the quantity b m such that, 

2~ b “” 2 < f - L f I! <2- b ”-'. 

K > z II oo 

Theorem 4.3.1: 

Let f be the original image and L r k 2 , m be the reconstruction. Let b m be the depth 
associated with m. Then the following holds: 

, b m C 

f e X k ^“ if and only if there exists a constant C such that ~ - a " “ ■ • 



Chapter 4 


125 


Proof; 


Let I* 


By theorem 3.4.2 we have. 

f-c. 

JU c 
11 




" f - mi 

2 b ”‘ : < 

|<- e, 

|j ^ mt 

l - j W j! ** 

mu 

-h,„- 

c 

.5- > 

c 

a - 



m m 


b C 

Conversely, let C be a constant such that — > a holds. 

m m 


=> nia ■ - b m < C 


[vf €X k+, ’°] 

[By definition 4.3.2] 



(• 2 >’» ! 

v i *• 

^ - nm 


[By theorem 3.4.3] 
[By definition 4.3.2] 


SC,2‘ 1 <C 3 (say) 

k.,(f.2 *)’ 


is unifonnly bounded sequence. :. f e X 


k+l,a 


From the above theorem we infer that, if for an image, the value of a is high then smaller 
value of m is sufficient for b-perfect reconstruction. This shows that for smooth images 
the storage requirement is less. This is in accordance with the experimental result (4.3.1). 
As the depth of test image decreases the image architecture starts deteriorating. We 
present a test image at various depths and we see that the deterioration is sharp as the 
number of bits per pixel (bpp) is reduced. 



Images at various depth 



Depth 5 bpp Depth 4 bpp 


Image 4.3.1 


In the next section we consider another important class of compression method, based on 


wavelets. 
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Chaptei 4 


4.4 Wavelet Basies: 

In this section, we present some wavelet preliminaries, wavelet based compression 
algorithm and experimental results. These are used for comparison with the results 

obtained through methods described in the previous section. 


Wavelets are elements generated from one single function \j/ by its dilation and 

translations. 


M ,ab (t) = |aT> 


1/2 


t - b 
V a 


oo 

where \p is the mother wavelet satisfying j vj/(x)dx = 0. In addition, v|/ has to satisfy 

“ oo 


r \Jf((0)| 

admissibility condition, which is J — t-t — do> < oo . High frequency wavelets correspond 

*00 * 

to a < 1 or narrow width, while low frequency wavelets have a > 1 or wider width. 


The main objective of the wavelet transform is to represent any arbitrary function as a 
superposition of wavelets. Such a superposition decomposes f into different scale levels, 
where each level is then further decomposed with a resolution adapted to the level. Let 
f e LjL-H ). A decomposition of f can be achieved by writing f as an integral over a and b 

ofV' h , 

W f(a,b) = |aj 1/2 |f(x)t|/^— Jdx. (4.4.1) 

”•00 

This refers to the continuous parameter wavelet transform, the dilation and translation 
parameters a and b run over 9? + . In most practical applications, discrete parameter 
wavelet transform is used. It is obtained by restricting the parameters a and b only to 
discrete values, a = a™ , b = nb 0 a™ , m, n e Z + , a 0 > 1, b 0 > 0 , fixed. 

oo 

W m , n (f) = C m ,„(f) =a- m/2 jf(t)i|/(ao m t -nb>. 

— 00 

Then the wavelet decomposition is given by, 


(4.4.2) 
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^ S ^ nun ( 0 V nun (4.4.3) 

nun 

V = a 0 m 2 i|/( a o m t - nb 0 ) . 

For a (l - 2 and b,i 1, there exist special choices of \|/ such that constitute an ortho 

normal basis and 

oc 

^ nun ( ^ ( l l / nun > = (x) \\f m Q (x) dx . (4.4.4) 

~00 

Different bases of this nature were constructed by Stromberg[70], Meyer[52], Battle[4] 
and Daubechies { 1 5], Mallat[49] invented a technique called Multi resolution analysis to 
the wavelet bases, in image analysis, and has given a fast computation algorithm. 

In Multi resolution analysis there are two functions, the mother wavelet vj/ and the scaling 
function <p. The dilated and translated versions of the scaling function 
f, u ,(x) - 2 m <j>(2 x n). For a fixed m, the cj) n%n are orthonormal. The space spanned 

by (j) m ,n is denoted by V,„ and the successive approximation spaces are of the form 
...Vj c V] c Vo c: V.j c V .2 ..., each with resolution 2 m . For each m, the span a 
space W m which is the orthogonal complement in V m .i of V m . The coefficients therefore 

m 1 

describe the information lost when going from an approximation of f with resolution 2 
to the coarser approximation with the resolution 2 m . These ideas are translated into the 
following algorithm for the computation of C m , n (f) (Daubechies [14], has more details): 

— X&2n-k ^m-l,k (-0 

k 

= (4.4.5) 

k 

Where g, = (-l)'h 1(I and h„ = 2 1/2 j<i>(x-n)cp(2x)dx . If the function f is given in 

sampled form, then one can take these samples for the highest order resolution 
approximation coefficients ao, n and the expressions for C m , n and a m>n describes a sub-band 
coding algorithm on these sampled values, with low pass filter h and high pass filter g. 
The reconstruction can be done using 

a „,- M (0 = X/ [^2n-i a m,n(f) + Szn-l^m.n (Oj • 


(4.4.6) 
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4.4.1 Compression Algorithm and storage bounds: 

Here we present the wavelet based compression strategy (R.A. DeVore et al., [17]). The 
two major steps iirvolvcd in compression using wavelet transform coding is to obtain the 
transform coefficients and then quantization of the transform coefficients. The quantized 

transform coefficents are stored instead of the pixel data. 


Step 1 : For the chosen wavelet basis, the coefficients C m , n f as in 4.4.3 are calculated 

using the pixel data. 

r iv 


Step 2: Choose a positive integer N, and a number 0 < a < A. Let x = 


a + — 

V 2 j 


The 


quantized coefficients C m n (f) , are chosen such that they satisfy the following: 

1 


(i) |C J , k {f)-C jJ .(f)| 

(ii) |C a | < C jlt , and 


(iii) 




N 


0 if C ;i <-. 

N 


j.r 


Step 3: The quantized coefficients C m n (f) ’s are stored using the entropy encoder based 
on Huffman coding (cf : 1.2.1). 


Using the bounds on the l p norm of the wavelet coefficients, in terms of the semi norm on 
Besov spaces, a bound on storage requirements can be obtained. The l p bound on wavelet 
coefficients for wavelets arising from a multi resolution analysis with the scaling function 
and the wavelet satisfying the following conditions are available in Wojtasazczyk [75]. 
Here we state it as a lemma. The conditions on scaling function and the wavelets are 

(i) |<Kx)| < C(l+ | x |)' A 

(ii) |<t>'(x)|<C(l+|x|)- A 

(iii) |\|/(x)| <C(l+|x|)" A 


for some A > 3 
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Lemma 4.4.1: I et a, 0 < a 
constant C such that. 


'A be given, and put 


x = 


l\ 

a + — 

2) 


Then there exists a 




where. If I 

7 i ii,a,T 


0 * ) 
sup t u o)(f;t) T 


1 < X < oo 

X = 00 


Theorem 4.4.1 (Storage Bound): 

Let r| be the number of non-zero coefficients C;j obtained after the application of 
Algorithm 4.4.1 for compression. Then r\ satisfies the following bound, 


ft 




Proof: 


Cj.k = (f,'| , j,k) anc * compressed coefficients C jk satisfy 
algorithm. We have, 


"j,k 


< Cj k as per the 


N 


i*k j,k i.k jx 


< c |fT 


iT,a,t 


[By lemma 4.4.1] 


Therefore we have. 


nscN|f|;„. 


4.4.2 Experiments: 

This is experimentally implemented on the test images. For experiments we use three 
wavelets, Daub4, Daubl2 and Daub20, developed by Daubechies[14]. The numbers 4,12 
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and 20 icpicsents the number of filter coefficients. We begin with an image which is 
coded at the rate of I byte per pixel. It has 256 gray levels. In these experiments, the 
quantized transform coefficients are stored. The quantization threshold determines the 
amount of compression. High frequency coefficients are normally made zero in the 
above algorithm. But the low frequency coefficients take more than a byte for storage. 
Since the nuinbei of coefficients, which are quantized, is sufficiently large, the net 
storage is lesser than the storage of MR image at 1 byte per pixel. - The results of 
algorithm 4.4. 1 implemented on img.60 (test image) are presented below. 


Img 60 using Daub4, Daub 12, Daub20 
Percentage compression with Relative Errors rll,rl2 


Daub 4 



65.71 67.29 68.26 69.08 69.62 70.06 70.45 


Graph 4.4.1 








Cha pter 4 


Daub 12 



Graph 4.4.2 


Daub 20 



Graph 4.4.3 
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COMPRESSION THROUGH WAVELETS (4.4.1) 
Images after decompression 



P.C.: 51.66 % B.R.: 3.87 bpp P.C.: 46.53 % B.R.: 4.277 bpp 


Image 4.4.1 






Compression through wavelets ( 4.4.1 ) 

Images after decompression 

Compression : 67 % Bitrate : 2.63 bpp 



Image 4.4.2 
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From the Images 4.4.1-4.4.2 and Graph 4.4.2 we note that ,h • 

information in the decompressed images if the amount of ” “ “ V ' SUal '° SS ^ 
than 56 % for the test imaees , m compresston done is not more 

-notation in the architect ^4^^^ " “ * 

.' 2 ' S bC " 1 Cr ,ha '; DaUl ' 4 md DaUb 20 fOT « i-ges considered here.TIpa“ 

images obtained after wavelet compression (Image 4 4 1 -4 4 21 and ■ , . 

v na S c -4.4.2) and images obtained after 

compressing through algorithm 4.2., (image 4 . 2 ., -4.2.5), we see tha, algorithm 4 2 
tmp emented using Poisson process is hotter sut.ed for MRI images. A comparison of the 

per otmance of standard eompre^ion techniques (Table 1.5.1) and that of algorithm 4.2.1 

implemented using Poisson process (421) for thp t*»et ; 

} r the test lma S es ^er consideration shows 
that the algorithm 4.2.1 performs better on the test images. 
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